Joint inversion method, device, storage medium and electronic equipment

Abstract

The application provides a joint inversion method and device, a storage medium and an electronic device. The method comprises the following steps: obtaining initial parameter data of a target area, and obtaining PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data comprises upper medium parameter data and lower medium parameter data; constructing a forward 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 a to-be-inverted parameter based on the target inversion function; and calculating target inversion results of each target parameter in a target parameter set based on the target inversion result of the to-be-inverted parameter, wherein the target parameter set comprises an anisotropy parameter. The embodiment of the application can improve the accuracy of seismic inversion.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oil and gas seismic exploration, and in particular to a joint inversion method and device, a storage medium and an electronic device. BACKGROUND

[0002] The change of pre-stack seismic inversion focused amplitude with offset or incidence angle is a key technology for describing reservoir elasticity, physical property and anisotropy; VTI medium (anisotropic medium) as a simplest anisotropic model usually refers to thin interbedded medium or uniform background with horizontal fractures, and plays an important role in the description of seismic wave propagation, and shale as a typical VTI medium model has obvious anisotropy characteristics, especially in the reflection and transmission response of seismic waves, that is, shale reservoirs usually exhibit strong anisotropy. However, the related technology usually focuses on the weak anisotropy assumption, and it is difficult to effectively predict anisotropy parameters from actual seismic data, resulting in low accuracy of seismic inversion. Based on this, how to improve the accuracy of seismic inversion currently does not have a good solution. SUMMARY

[0003] Therefore, the embodiments of the present application provide a joint inversion method and device, a storage medium and an electronic device to solve the problem of low accuracy of seismic inversion in related technologies. That is, the embodiments of the present application can realize a PP wave and PS wave joint inversion method through PP wave seismic data and PS wave seismic data, to predict the anisotropy parameters of VTI medium, so as to obtain a target inversion result of anisotropy parameters with high accuracy, and effectively improve the accuracy of seismic inversion.

[0004] According to an aspect of the embodiments of the present application, a joint inversion method is provided, and the method comprises:

[0005] Obtaining initial parameter data of a target area, and obtaining PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data comprises upper medium parameter data and lower medium parameter data;

[0006] 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;

[0007] Based on the target inversion function, a target inversion result of a to-be-inverted parameter is determined;

[0008] Based on the target inversion result of the to-be-inverted parameter, a target inversion result of each target parameter in a target parameter set is calculated, and the target parameter set comprises anisotropy parameters.

[0009] According to another aspect of the embodiments of the present application, there is provided a joint inversion device, the device comprising:

[0010] an acquisition unit configured to acquire initial parameter data of a target region, and to acquire PP wave seismic data and PS wave seismic data of the target region, the initial parameter data comprising upper medium parameter data and lower medium parameter data;

[0011] 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;

[0012] the processing unit is further configured to determine a target inversion result of a parameter to be inverted based on the target inversion function;

[0013] the processing unit 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, the target parameter set comprising anisotropy parameters.

[0014] According to another aspect of the embodiments of the present application, there is provided an electronic device comprising a processor and a memory storing a program, wherein the program comprises instructions which, when executed by the processor, cause the processor to perform the above-mentioned method.

[0015] According to another aspect of the embodiments of the present application, there is provided a non-transitory computer-readable storage medium storing computer instructions for causing a computer to perform the above-mentioned method.

[0016] According to another aspect of the embodiments of the present application, there is provided a computer program product comprising a computer program which, when executed by a processor, is configured to cause a computer to perform the above-mentioned method.

[0017] The embodiments of the present application can, after acquiring initial parameter data of a target region and acquiring PP wave seismic data and PS wave seismic data of the target region, 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, the target inversion result of a parameter to be inverted can be determined based on the target inversion function, and the target inversion result of each target parameter in a target parameter set can be calculated based on the target inversion result of the parameter to be inverted, the target parameter set comprising anisotropy parameters. It can be seen that the embodiments of the present application can realize a PP wave and PS wave joint inversion method through the PP wave seismic data and the PS wave seismic data, so as to predict anisotropy parameters of a VTI medium, thereby obtaining a target inversion result of anisotropy parameters with high accuracy, and effectively improving the accuracy of seismic inversion. BRIEF DESCRIPTION OF DRAWINGS

[0018] In the following description of the example embodiments in conjunction with the drawings, more details, features and advantages of the present application are disclosed, in which:

[0019] Figure 1 A flowchart of a joint inversion method according to an example embodiment of the present application is shown;

[0020] Figure 2 A flowchart of another joint inversion method according to an example embodiment of the present application is shown;

[0021] Figure 3 A schematic diagram of a PP-wave reflection coefficient according to an example embodiment of the present application is shown;

[0022] Figure 4 A schematic diagram of a PS-wave reflection coefficient according to an example embodiment of the present application is shown;

[0023] Figure 5 A schematic diagram of an inversion result according to an example embodiment of the present application is shown;

[0024] Figure 6 A schematic diagram of another inversion result according to an example embodiment of the present application is shown;

[0025] Figure 7 A schematic diagram of an inversion result of an anisotropic parameter according to an example embodiment of the present application is shown;

[0026] Figure 8 A schematic diagram of another inversion result of an anisotropic parameter according to an example embodiment of the present application is shown;

[0027] Figure 9 A schematic diagram of yet another inversion result of an anisotropic parameter according to an example embodiment of the present application is shown;

[0028] Figure 10 A schematic diagram of still another inversion result of an anisotropic parameter according to an example embodiment of the present application is shown;

[0029] Figure 11 A schematic block diagram of a joint inversion apparatus according to an example embodiment of the present application is shown;

[0030] Figure 12 A structural block diagram of an example electronic device that can be used to implement embodiments of the present application is shown. DETAILED DESCRIPTION

[0031] Embodiments of the present application will be described below in greater detail with reference to the accompanying drawings. While certain embodiments of the present application are shown in the drawings, it is understood that the present application can be embodied in various forms and should not be interpreted in limitation to the embodiments set forth herein, but rather by way of example. It is to be understood that the drawings and embodiments are only for illustrative purposes and are not intended to limit the scope of the present application.

[0032] It should be understood that each of the steps in the method embodiments of the present application can be performed in a different order and / or in parallel. In addition, the method embodiments can include additional steps and / or omit performing the steps shown. The scope of the present application is not limited in this respect.

[0033] The term "comprising" and variations thereof as used herein are used inclusively, i.e., "comprising, but not limited to." The term "based on" is "based, at least in part, on." The term "one embodiment" means "at least one embodiment"; the term "another embodiment" means "at least one additional embodiment"; the term "some embodiments" means "at least some embodiments." Related terms are defined in the description that follows. It should be noted that reference herein to "first", "second", etc. concepts is merely used to differentiate different apparatuses, modules, or units, and is not intended to limit the order or interdependence of the functions performed by these apparatuses, modules, or units.

[0034] It should be noted that the terms "a" and "an" and "the" and similar referents used in the context of describing the application are to be construed to be inclusive, while also including the singular, unless the context dictates otherwise. It should be noted that the terms "first", "second", etc. are used herein merely to differentiate one element from another, and are not intended to limit the order or interdependence of the functions performed by these elements.

[0035] The names of the messages or information exchanged between the various apparatuses in the embodiments of the present application are used merely for illustrative purposes, and are not intended to limit the scope of these messages or information.

[0036] It should be noted that the execution subject of the joint inversion method provided by the embodiments of the present application can be one or more electronic devices, and the present application does not limit this; wherein the electronic device can be a terminal (i.e. client) or a server, so 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 by the embodiments of the present application can be executed by the terminal and the server together. Correspondingly, the terminal mentioned here can include but not limited to: smart phone, notebook computer, desktop computer, smart voice interaction device, etc. The server mentioned here can be a standalone physical server, or a server cluster or distributed system composed of multiple physical servers, or a cloud server providing cloud service, cloud database, cloud computing (cloud computing), cloud function, cloud storage, network service, cloud communication, middleware service, domain name service, security service, CDN (Content Delivery Network, content distribution network), and basic cloud computing services such as big data and artificial intelligence platform, etc.

[0037] Based on the above description, the embodiments of the present application propose a joint inversion method, which can be executed by the above-mentioned electronic device (terminal or server); or the joint inversion method can be executed by the terminal and the server together. In order to facilitate description, the joint inversion method executed by the electronic device is taken as an example for description hereinafter; for example, as shown in the figure, the joint inversion method can include the following steps S101-S104: Figure 1

[0038] S101, 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 including upper medium parameter data and lower medium parameter data.

[0039] Optionally, the target region can be any region, and the embodiments of the present application do not limit this.

[0040] ​In the embodiments of the present application, the initial parameter data can include initial values of each target parameter in the target parameter set; optionally, the target parameter set can include, but is not limited to, at least one of the following: a transverse wave velocity (i.e., qS wave (also referred to as S wave) velocity, qS wave is a transverse wave), a longitudinal wave velocity (i.e., qP wave (also referred to as P wave) velocity, qP wave is a longitudinal wave), a density, and an anisotropy parameter, etc., and the embodiments of the present application do not make any limitation in this regard. Optionally, the anisotropy parameter can 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., and the embodiments of the present application do not make any limitation in this regard. Optionally, the initial value of one target parameter can include an upper initial value of the corresponding target parameter in the upper medium parameter data and / or a lower initial value of the corresponding target parameter in the lower medium parameter data. Optionally, the initial parameter data can further include an angle between a ray and a normal of an interface, etc., in which case the following vertical slowness, etc., can be determined based on the initial parameter data; or, the initial parameter data can further include a longitudinal wave vertical slowness, a transverse wave vertical slowness, and a horizontal slowness, etc., and the embodiments of the present application do not make any limitation in this regard.

[0041] In the embodiments of the present application, the PP wave can refer to a P wave as incident and a P wave as reflected (PP wave can also be referred to as reflected longitudinal wave), and the PS wave can refer to a P wave as incident and an S wave as reflected (PS wave can also be referred to as converted reflected wave).

[0042] Optionally, the initial parameter data, the PP wave seismic data, and the PS wave seismic data can be data in a logging model, i.e., can be generated in model testing, or can be data in actual seismic data, etc., and the embodiments of the present application do not make any limitation in this regard.

[0043] Optionally, the initial parameter data, the PP wave seismic data, and the PS wave seismic data of the target region can be acquired in the following ways, including but not limited to the following:

[0044] The first acquisition method: the initial parameter data, the PP wave seismic data, and the PS wave seismic data of each region in a plurality of regions can be stored in the storage space of the electronic device, and the plurality of regions include the target region, in which case the initial parameter data, the PP wave seismic data, and the PS wave seismic data of the target region can be acquired from the storage space.

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

[0046] The third obtaining manner: the electronic device can display a 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 of the target region, the PP wave seismic data and the PS wave seismic data through the parameter input operation, in this case, the electronic device can respond to the parameter input operation and obtain the initial parameter data of the target region, the PP wave seismic data and the PS wave seismic data, and the like.

[0047] In S102, a forward operator is constructed based on the initial parameter data, and a target inversion function is constructed based on the forward operator, the PP wave seismic data and the PS wave seismic data.

[0048] Optionally, the forward operator can include a first forward operator and a second forward operator; optionally, the first forward operator can also be referred to as a PP wave forward operator, and the second forward operator can also be referred to as a PS wave forward operator.

[0049] In the embodiment of the application, when the forward operator is constructed based on the initial parameter data, the electronic device can determine a stiffness matrix based on the initial parameter data, the stiffness matrix can include a plurality of stiffness elements (i.e., one stiffness matrix can include element values of each stiffness element in the plurality of stiffness elements; optionally, the element value of one stiffness element can be an initial element value of the corresponding stiffness element), the initial element value of one 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 the initial values of each target parameter in the target parameter set, and the initial value of one target parameter includes an upper initial value of the corresponding target parameter in the upper medium parameter data and / or a lower initial value of the corresponding target parameter in the lower medium parameter data. Then, a plurality of PP wave parameter perturbation coefficients can be calculated based on the stiffness matrix, and a plurality of PS wave parameter perturbation coefficients can be calculated based on the stiffness matrix; wherein the plurality of PP wave parameter perturbation coefficients include 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 PS wave parameter perturbation coefficients corresponding to each stiffness element; optionally, the plurality of PP wave parameter perturbation coefficients can also include a PP wave density parameter perturbation coefficient (i.e., the fifth PP wave parameter perturbation coefficient described below), and the plurality of PS wave parameter perturbation coefficients can also include a PS wave density parameter perturbation coefficient (i.e., the fifth PS wave parameter perturbation coefficient described below). Based on this, the electronic device can construct the first forward operator based on the plurality of PP wave parameter perturbation coefficients, and construct the second forward operator based on the plurality of PS wave parameter perturbation coefficients, to realize the construction of the forward operator.

[0050] It should be noted that the elastic impedance tensor can accurately represent the reflection and transmission coefficients at any incident angle and any elastic contrast interface, and the elastic impedance tensor of the VTI medium can be represented as shown in formula 1.1:

[0051] Formula 1.1

[0052] wherein, E represents an elastic impedance tensor (such as E P may represent a longitudinal wave elastic impedance tensor, E S may represent a transverse wave elastic impedance tensor), the subscripts P and S may respectively represent qP wave and qS wave, and p and q may respectively represent horizontal slowness and vertical slowness (such as q P may represent longitudinal wave vertical slowness, q S may represent transverse wave vertical slowness), C may represent a stiffness matrix, C 55 may represent stiffness elements of the stiffness matrix.

[0053] In an embodiment of the present application, the expression of each stiffness element in the stiffness matrix may be as shown in Formula 1.2, that is, the electronic device may adopt Formula 1.2 to calculate the stiffness elements in the stiffness matrix, that is, to calculate the element value of each stiffness element in the stiffness elements.

[0054] Formula 1.2

[0055] wherein, V P may represent longitudinal wave velocity, V S may represent transverse wave velocity, ρ may represent density, δ may represent a longitudinal anisotropy parameter (i.e., a first anisotropy parameter), ε may represent a transverse anisotropy parameter (i.e., a second anisotropy parameter), and γ may represent a y-direction anisotropy parameter (i.e., a third anisotropy parameter); f = 1 - C 55 / C 33 .

[0056] Based on this, the electronic device may 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 may be the parameter data of the upper medium, and the lower medium parameter data may 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 longitudinal wave velocity, transverse wave velocity, and 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 an anisotropy parameter (i.e., a medium is an isotropic medium), the value of the anisotropy parameter of the medium may be 0.

[0057] Based on this, the electronic device can determine the upper layer initial value (i.e., the upper layer medium initial value) and the lower layer initial value (i.e., the lower layer medium initial value) of each target parameter based on the initial parameter data, thereby realizing the determination of the initial value of each target parameter, and performing mean value operation on the upper layer initial value and the lower layer 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 value between the upper layer initial value and the lower layer initial value of the corresponding target parameter, at this time the parameter value of a target parameter can also be called the initial parameter value), and using the parameter value of each target parameter to determine the initial element value of each stiffness element in the plurality of stiffness elements (i.e., at this time the initial parameter value of each target parameter can be used to determine the initial element value of each stiffness element), so as to realize the determination of the stiffness matrix. For example, 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, that is, 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 can include but are not limited to: a first stiffness element (i.e., C 11 ), a second stiffness element (i.e., C 13 ), a third stiffness element (i.e., C 33 ), and a fourth stiffness element (i.e., C 55 ), etc., which are not limited by the embodiments of the present application.

[0058] Further, when calculating the plurality of PP wave parameter perturbation coefficients based on the stiffness matrix, the electronic device can determine a plurality of intermediate coefficients based on the stiffness matrix; and calculate the PP wave diagonal parameter based on the plurality of intermediate coefficients, the longitudinal wave polarization direction coefficient, and the longitudinal wave vertical slowness; based on this, the plurality of PP wave parameter perturbation coefficients can be calculated based on the PP wave diagonal parameter, the stiffness matrix, and the plurality of intermediate coefficients. Correspondingly, when calculating the plurality of PS wave parameter perturbation coefficients based on the stiffness matrix, the PS wave diagonal parameter can be calculated based on the plurality of intermediate coefficients, the transverse wave polarization direction coefficient, and the transverse wave vertical slowness; and the plurality of PS wave parameter perturbation coefficients can be calculated based on the PS wave diagonal parameter, the stiffness matrix, and the plurality of intermediate coefficients.

[0059] It should be noted that in an anisotropic medium, the phase velocity and the polarization direction change with the propagation direction. The polarization direction of the VTI medium can be represented as formula 1.3:

[0060] Formula 1.3

[0061] Wherein, U can represent the polarization direction (such as U P can represent the longitudinal wave polarization direction, U S can represent the transverse wave polarization direction), 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 Pmay represent a coefficient of polarization direction of a longitudinal wave, G S may represent a coefficient of polarization direction of a transverse wave); it should be understood that the polarization direction x components of the qP wave and the qS wave are both p, i.e., the same. Wherein, the vertical slowness of the longitudinal wave can be calculated by the longitudinal wave velocity (such as calculating the vertical slowness of the longitudinal wave by using the parameter value of the longitudinal wave velocity (such as taking 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 as the vertical slowness of the longitudinal wave) to realize the initialization of the vertical slowness of the longitudinal wave), and the vertical slowness of the transverse wave can be calculated by the transverse wave velocity (such as calculating the vertical slowness of the transverse wave by using the parameter value of the transverse wave velocity (such as taking 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 as the vertical slowness of the transverse wave) to realize the initialization of the vertical slowness of the transverse wave); 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 realize the initialization of the horizontal slowness, and so on. Based on this, the embodiments of the present application can initialize the stiffness matrix, the vertical slowness of the longitudinal wave, the vertical slowness of the transverse wave, and the horizontal slowness, and so on, by the parameter values of the various target parameters (i.e., by the initial parameter data). Accordingly, the electronic device can calculate the coefficient of polarization direction of the longitudinal wave based on the parameter value of the density, the horizontal slowness, the vertical slowness of the longitudinal wave, and the stiffness matrix (i.e., the element value of at least one stiffness element in the plurality of stiffness elements), and can calculate the coefficient of polarization direction of the transverse wave based on the parameter value of the density, the horizontal slowness, the vertical slowness of the transverse wave, and the stiffness matrix; that is, the electronic device can calculate the coefficients of polarization direction of the longitudinal wave and the transverse wave by using the calculation formula of the polarization direction coefficient in formula 1.3.

[0062] Based on this, according to the elastic impedance tensor and the polarization direction of the VTI medium, the propagation matrix can be represented as shown in formula 1.4:

[0063] Formula 1.4

[0064] Under the boundary conditions of displacement continuity and stress continuity, the exact reflection coefficient can be represented by the polarization direction and the propagation matrix as shown in formula 1.5:

[0065] Formula 1.5

[0066] Wherein, R=[R PP R SP ;R PS R SSis 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) can represent the upper medium and the lower medium respectively, and the subscripts "+" and "-" can represent the down-going wave and the up-going wave respectively. Considering that the vertical slowness of the up-going wave and the down-going wave is opposite, according to the formula 1.3 and the formula 1.4, the relationship between the polarization direction of the up-going wave and the down-going wave and the propagation matrix of the up-going wave and the down-going wave can be shown as the formula 1.6:

[0067] Formula 1.6

[0068] Based on this, the formula 1.6 is substituted into the formula 1.5, and the accurate reflection coefficient can be simplified as shown in the formula 1.7:

[0069] Formula 1.7

[0070] Wherein, the matrix A and B can be represented as shown in the formula 1.8:

[0071] Formula 1.8

[0072] Wherein, the superscript (m) (m=1,2) represents the upper medium and the lower medium, and I=[1 0;0 1] is the unit matrix. Based on this, the embodiment of the application can realize the accurate reflection coefficient of the up-going wave and the down-going wave by the 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 expression of the matrix A and B can be shown as the formula 1.9:

[0073] Formula 1.9

[0074] Based on this, the embodiment of the application constructs the accurate expression of the VTI medium reflection coefficient represented by two matrices A and B (the scale is 2*2). Correspondingly, according to the specific expression of the matrix A and B (i.e. formula 1.9), the formula 1.10 can be obtained:

[0075] Formula 1.10

[0076] Wherein, Λ=diag([Λ P ,Λ S ]) is a diagonal matrix, Λ P (i.e. PP wave diagonal parameter) can be the result of substituting the first column of the matrix A and the first column of the matrix B into the formula 1.10 for calculation, and Λ S(i.e., the PS wave diagonal parameter) can be substituted into formula 1.10 for the second column of matrix A and the second column of matrix B to calculate the result. According to formula 1.10, the inverse operation of matrix A and B can be converted into a transpose operation, and under the weak contrast assumption condition (i.e., the weak contrast assumption condition 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.), 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:

[0077] Formula 1.11

[0078] It should be understood that petrophysical experiments show 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 application can use the weak contrast assumption, that is, the elastic parameters and anisotropy parameters of the upper and lower media are small, which can improve the accuracy of the reflection coefficient at large angles in the strong anisotropy medium, that is, the conventional weak anisotropy assumption eliminates the high-order terms of the anisotropy parameters, or the difference is not much, such as only the first-order term is retained, which causes the accuracy to be damaged, especially at large angles. The strong anisotropy parameter assumption mentioned in the embodiment of the present application can effectively improve the accuracy, especially the accuracy of the reflection coefficient at large angles in the strong anisotropy medium. Moreover, formula 1.11 gives a general approximation formula of the reflection coefficient that can be applied to strong anisotropy VTI medium, which can be simply represented by two 2x2 matrices and transpose operation, thereby effectively simplifying the calculation and improving the operation efficiency.

[0079] In addition, the conventional reflection coefficient approximation equation, which is usually represented by elastic parameters and anisotropy parameter perturbation, has a problem of great difference in contribution between parameters, which seriously affects the stability of pre-stack seismic inversion, especially the anisotropy parameter; and the general reflection coefficient approximation formula (i.e., formula 1.11) derived by the embodiment of the application can be represented by the stiffness matrix, which overcomes the problem of great difference in contribution between parameters; that is, the conventional method uses five parameters for calculation, and the anisotropy is not very large even if the anisotropy parameter is large; and the wave velocity (such as the longitudinal wave velocity and the transverse wave velocity) is large, and the contribution to the reflection coefficient is much larger; the parameters with large contribution are better inverted, and the parameters with small contribution are not easy to invert, which leads to that the inversion of the anisotropy parameter is difficult and the precision is low; and the embodiment of the application binds the anisotropy parameter and the wave velocity together through the stiffness matrix, which can effectively reduce the difference, thereby facilitating seismic inversion, and effectively improving the inversion precision of the anisotropy parameter.

[0080] Further, the embodiment of the application can derive the weak-contrast reflection coefficient approximation of PP wave and PS wave of VTI medium based on the general reflection coefficient approximation formula. Based on this, formula 1.9 can be substituted into formula 1.11 to obtain the reflection coefficient equation of PP wave and PS wave of VTI medium; wherein the PP wave reflection coefficient can be represented as formula 1.12:

[0081] Formula 1.12

[0082] Wherein, 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 , ; optionally, may be the stiffness element difference between the upper and lower media (Y can be C 11 , C 13 , C 33 , and C 55), such as the upper-layer element value of each stiffness element can be determined by using the upper-layer initial value of each target parameter, the lower-layer element value of each stiffness element can be determined by using the lower-layer initial value of each target parameter, and the difference between the lower-layer element value and the upper-layer element value of each stiffness element can be taken as the stiffness element difference of the corresponding stiffness element (i.e., the stiffness element difference of the corresponding stiffness element between the upper-layer medium and the lower-layer medium), and correspondingly, 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, and the like; optionally, L can be the ratio of the transverse wave velocity to the longitudinal wave velocity. Based on this, the electronic device can determine a plurality of intermediate coefficients based on the stiffness matrix, which can include Γ 11 (also referred to as the first intermediate coefficient), Γ 13 (also referred to as the second intermediate coefficient), Γ 33 (also referred to as the third intermediate coefficient), Γ 55 (also referred to as the fourth intermediate coefficient), and Γ ρ (also referred to as the density intermediate coefficient). Wherein, Λ P may be as shown in formula 1.13:

[0083] Formula 1.13

[0084] Based on this, the electronic device can calculate the PP wave diagonal parameter based on the plurality of intermediate coefficients, the longitudinal wave polarization direction coefficient, and the longitudinal wave vertical slowness through formula 1.13. Correspondingly, ξ P (also referred to as the first reflection coefficient calculation parameter) can be as shown in formula 1.14:

[0085] Formula 1.14

[0086] Correspondingly, Q ψ (also referred to as the second reflection coefficient calculation parameter) can be as shown in formula 1.15:

[0087] Formula 1.15

[0088] In addition, the VTI medium PS wave reflection coefficient can be expressed as shown in formula 1.16:

[0089] Formula 1.16

[0090] Wherein, Λ S (i.e., the PS wave diagonal parameter) can be as shown in formula 1.17:

[0091] Formula 1.17

[0092] Based on this, the electronic device can calculate the PS wave diagonal parameter based on the plurality of intermediate coefficients, the transverse wave polarization direction coefficient, and the transverse wave vertical slowness through formula 1.17. Correspondingly, The third reflection coefficient calculation parameter (also referred to as the third reflection coefficient calculation parameter) can be as shown in formula 1.18:

[0093] Formula 1.18

[0094] Correspondingly, ξ S The fourth reflection coefficient calculation parameter (also referred to as the fourth reflection coefficient calculation parameter) can be as shown in formula 1.19:

[0095] Formula 1.19

[0096] It is worth noting that there is a ratio of horizontal and vertical slowness in the coefficient term, and the horizontal and vertical slowness are the ratio 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 is irrelevant to the phase velocity, and only relevant to the anisotropy parameter and the transverse-longitudinal wave velocity ratio, that is, the present application implements the derivation of a PP wave and PS wave linear reflection coefficient.

[0097] Further, considering that the stiffness matrix is used to calculate the anisotropy parameter, the embodiment of the present application can construct a prestack seismic inversion target functional (that is, a target inversion function) constrained by the relationship between the stiffness matrix and the anisotropy parameter. It should be understood that using conventional prestack seismic inversion methods, the stability of directly predicting the anisotropy parameter from seismic data is insufficient (because the contribution of the anisotropy parameter is small); and the reflection coefficient equation derived by the embodiment of the present application can be represented by the stiffness matrix, which overcomes the problem that the contribution of each parameter is different in the conventional reflection coefficient approximation, making it possible to predict the stiffness matrix from seismic data and indirectly predict the anisotropy parameter stably.

[0098] In the embodiment of the present application, the PP wave and PS wave reflection coefficient approximation formula (that is, formula 1.12 and formula 1.16) can be represented by the perturbation terms of the stiffness coefficients (that is, the stiffness elements) C 11 , C 13 , C 33 , C 55 and ρ, and these perturbation terms are written in the form of logarithmic difference as shown in formula 1.20:

[0099] Formula 1.20

[0100] Based on this, the forward model of the PP wave and the PS wave can be represented as shown in formula 1.21:

[0101] Formula 1.21

[0102] Where dP =[d P (θ P1 ),d P (θ P2 ),…,d P (θ Ph )] T may represent PP wave seismic data composed of h angle gathers with incident angles θ P1 , θ P2 , …, θ Ph , d S =[d S (θ S1 ),d P (θ S2 ),…,d P (θ Sh )] T may represent PS wave seismic data composed of h angle gathers with incident angles θ S1 , θ S2 , …, θ Sh , X=[X 11 ,X 13 ,X 33 ,X 55 ,ρ] T may be parameters to be inverted, the first forward operator G1 and the second forward operator G2 in the forward operator may be shown in formula 1.22:

[0103] Formula 1.22

[0104] wherein W P and W S may respectively represent wavelets (such as amplitude phase wavelets) extracted from the PP wave seismic data and the PS wave seismic data, D is a first order partial differential operator, v and v’ may be respectively a PP wave parameter perturbation coefficient and a PS wave parameter perturbation coefficient; optionally, the electronic device may use formula 1.23 to calculate the parameter perturbation coefficients (i.e. multiple PP wave parameter perturbation coefficients (which may include a first PP wave parameter perturbation coefficient v 11 , a second PP wave parameter perturbation coefficient v 13 , a third PP wave parameter perturbation coefficient v 33 , a fourth PP wave parameter perturbation coefficient v 55 and a fifth PP wave parameter perturbation coefficient v ρ ) and multiple PS wave parameter perturbation coefficients (which may include a first PS wave parameter perturbation coefficient , a second PS wave parameter perturbation coefficient , a third PS wave parameter perturbation coefficient , a fourth PS wave parameter perturbation coefficient and a fifth PS wave parameter perturbation coefficient )):

[0105] Formula 1.23

[0106] Based on this, the electronic device can calculate the plurality of PP wave parameter perturbation coefficients based on the PP wave diagonal parameter, the stiffness matrix and the plurality of intermediate coefficients through Formula 1.23; specifically, the electronic device can calculate the plurality of second reflection coefficient calculation parameters (i.e., Q ψ ) based on the parameter values of the stiffness matrix and the density through Formula 1.15, so as to calculate the plurality of PP wave parameter perturbation coefficients by using the PP wave diagonal parameter, the plurality of intermediate coefficients and the plurality of second reflection coefficient calculation parameters. Correspondingly, the electronic device can calculate the plurality of PS wave parameter perturbation coefficients based on the PS wave diagonal parameter, the stiffness matrix and the plurality of intermediate coefficients through Formula 1.23; in this case, after determining the plurality of intermediate coefficients based on the stiffness matrix, the electronic device can calculate the third reflection coefficient calculation parameter and the fourth reflection coefficient calculation parameter based on the plurality of intermediate coefficients and the transverse wave polarization direction coefficient respectively through Formula 1.18 and Formula 1.19, so as to calculate the plurality of PS wave parameter perturbation coefficients by using the PS wave diagonal parameter, the plurality of intermediate coefficients, the plurality of second reflection coefficient calculation parameters, the third reflection coefficient calculation parameter and the fourth reflection coefficient calculation parameter.

[0107] In the embodiment of the present application, when the first forward operator is constructed based on the plurality of PP wave parameter perturbation coefficients, the first wavelet (i.e., W P ) can be extracted from the PP wave seismic data (which can include h angle seismic data, h being a positive integer), and the first forward operator is constructed based on the first wavelet and the plurality of PP wave parameter perturbation coefficients, such as the calculation method of G1 in Formula 1.22. Optionally, when the first wavelet is extracted, each angle seismic data in the PP wave seismic data can be stacked to obtain PP wave stacked seismic data, so as to extract the first wavelet from the PP wave stacked seismic data; or the wavelet corresponding to each angle seismic data can be extracted from each angle seismic data in the PP wave seismic data respectively, and the mean operation (i.e., the mean operation is performed on the value of each sampling point) is performed on the wavelet corresponding to each angle seismic data to obtain the first wavelet, so as to extract the first wavelet from the PP wave seismic data, and the like; the present application is not limited thereto.

[0108] Correspondingly, when the second forward operator is constructed based on the plurality of PS wave parameter perturbation coefficients, the second wavelet (i.e., W S), and a second forward operator is constructed based on the second wavelet and the plurality of PS wave parameter perturbation coefficients, such as the second forward operator is constructed by the calculation manner of G2 in formula 1.22. It should be noted that the implementation manner of extracting the second wavelet can be the same as the implementation manner of extracting the first wavelet, and the embodiment of the present application will not be repeated here.

[0109] It should be understood that the coefficient terms in formula 1.23 are functions of anisotropy parameters and P-S wave velocity ratio, so the forward model is linear. In order to overcome the ill-posedness of calculating anisotropy parameters by using the predicted stiffness matrix in pre-stack seismic inversion, the nonlinear relationship between the stiffness matrix and the anisotropy parameters can be shown as formula 1.24:

[0110] Formula 1.24

[0111] In formula 1.24, the function g can represent the function of calculating anisotropy parameters by using the stiffness matrix, that is, each term on the right side of formula 1.24 can be used to calculate V P , V S , ρ, δ, ε in turn. Based on this, according to the PP wave and PS wave forward model, a target functional (i.e. target inversion function) constrained by the relationship between the stiffness matrix and the anisotropy parameters is constructed as shown in formula 1.25, that is, the electronic device can use formula 1.25 to construct the target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data:

[0112] Formula 1.25

[0113] In formula 1.25, 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, and the third regularization term can be the stiffness matrix and anisotropy parameter relationship constraint term. The regularization parameter α is a constraint weight, and the parameter α has less influence on the stiffness matrix inversion result, but has greater influence on the final anisotropy parameters. The greater the α is, the stronger the stability of the anisotropy parameter prediction is. Optionally, the regularization parameter can be set according to experience, can be set according to actual demand, can be obtained according to experimental debugging, etc. The embodiment of the present application does not limit this.

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

[0115] In the embodiment of the present application, the electronic device can perform inversion on the to-be-inverted parameter through the target inversion function, and then the target inversion result of the to-be-inverted parameter can be obtained.

[0116] In the embodiment of the present application, the electronic device can perform inversion on the to-be-inverted parameter through the target inversion function, and then the target inversion result of the to-be-inverted parameter can be obtained.

[0117] In the embodiment of the present application, the electronic device can perform inversion on the to-be-inverted parameter through the target inversion function, and then the target inversion result of the to-be-inverted parameter can be obtained.

[0118] In the embodiment of the present application, the electronic device can determine the calculation mode corresponding to each target parameter in the target parameter set based on the nonlinear relationship (which can be shown in formula 1.24) between the stiffness matrix and the target parameter set; and calculate the target inversion result of each target parameter according to the calculation mode corresponding to each target parameter and the target inversion result of the to-be-inverted parameter. Optionally, at this time, the electronic device can determine the calculation mode corresponding to the P wave velocity, the S wave velocity, the density, the first anisotropic parameter and the second anisotropic parameter respectively, so as to calculate the target inversion result of the P wave velocity, the S wave velocity, the density, the first anisotropic parameter and the second anisotropic parameter respectively; that is, the target inversion result of the to-be-inverted parameter X can be substituted into formula 1.24, so as to calculate the target inversion result of each target parameter according to the calculation mode corresponding to each target parameter and the target inversion result of the to-be-inverted parameter. Exemplarily, the target inversion result of the P wave velocity can be equal to sqrt(exp(X 33 ) / exp(X ρ )), the target inversion result of the S wave velocity can be equal to sqrt(exp(X 55 ) / exp(X ρ )), the target inversion result of the density can be equal to exp(X ρ ), and so on; that is, exp(X 11 ) can be taken as C 11 , exp(X 13 ) can be taken as C 13 , and exp(X 33 ) can be taken as C 33, exp(X 55 ) as C 55 , and exp(X ρ ) as p, substitute into formula 1.24, so as to calculate the target inversion result of each target parameter respectively, that is, the target inversion result of the longitudinal wave velocity, the target inversion result of the transverse wave velocity, the target inversion result of the density, the target inversion result of the first anisotropy parameter and the target inversion result of the second anisotropy parameter are sequentially calculated.

[0119] After the initial parameter data of the target region is acquired, and the PP wave seismic data and the PS wave seismic data of the target region are acquired, the embodiment of the present application can construct the forward operator based on the initial parameter data, and construct the target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data. Further, the target inversion result of the to-be-inverted parameter can be determined based on the target inversion function, and the target inversion result of each target parameter in the target parameter set, which includes the anisotropy parameter, can be calculated based on the target inversion result of the to-be-inverted parameter. It can be seen that the embodiment of the present application can realize the PP wave and PS wave joint inversion method through the PP wave seismic data and the PS wave seismic data, so as to realize the prediction of the anisotropy parameter of the VTI medium, thereby obtaining the target inversion result of the anisotropy parameter with higher accuracy, and effectively improving the accuracy of the seismic inversion.

[0120] Based on the above description, the embodiment of the present application further proposes a more specific joint inversion method. Correspondingly, 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. In order to facilitate description, the joint inversion method executed by the electronic device is taken as an example for description hereinafter; please refer to Figure 2 , the joint inversion method can include the following steps S201-S208:

[0121] S201, acquiring the initial parameter data of the target region, and acquiring the PP wave seismic data and the PS wave seismic data of the target region, the initial parameter data including the upper medium parameter data and the lower medium parameter data.

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

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

[0124] Optionally, in the 0th iteration (where k represents the iteration number, k=0, k is a positive integer), the electronic device can determine the current parameter data of the parameter to be inverted (i.e., the value of X in the kth iteration) by using the initial parameter data; that is, the current parameter data of the parameter to be inverted can be determined by using the parameter values of the target parameters determined by the initial parameter data to initialize the parameter to be inverted; for example, the current parameter data of the parameter to be inverted can be determined based on the calculation method corresponding to each target parameter and the parameter values of the target parameters, and the like.

[0125] Optionally, when calculating the partial derivative of the target inversion function based on the current parameter data of the parameter to be inverted, the electronic device can calculate the first-order partial derivative and the second-order partial derivative of the target inversion function based on the current parameter data of the parameter to be inverted to calculate the partial derivative of the target inversion function. Optionally, the electronic device can calculate the first-order partial derivative and the second-order partial derivative of the target inversion function by using formula 2.1:

[0126] Formula 2.1

[0127] wherein, and may represent the first-order partial derivative and the second-order partial derivative, respectively, and J and H can be the Jacobian matrix and the Hessian matrix of the function g, which can be specifically as shown in formula 2.2:

[0128] Formula 2.2

[0129] Based on this, the electronic device can calculate the Jacobian matrix and the Hessian matrix based on the current parameter data of the parameter to be inverted to calculate the first-order partial derivative and the second-order partial derivative of the target inversion function through the Jacobian matrix and the Hessian matrix.

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

[0131] In the embodiments of the present application, the iterative Newton method can be used for solving; according to the iterative Newton method, the update 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 update the parameter to be inverted by using formula 2.3 to obtain the updated parameter data of the parameter to be inverted:

[0132] Formula 2.3

[0133] wherein, k is the iteration number, and the parameter μ k is the iteration step size in the kth iteration, and X k represents the value of the parameter to be inverted in the kth iteration (i.e., the current parameter data of the parameter to be inverted), and X k+1The updated inversion parameter data can represent the to-be-inverted parameter; optionally, the iteration step length at any two iteration times can be the same or different, and embodiments of the present application do not limit this; for example, the iteration step length at each iteration time can be set to a specified step length according to experience or actual needs, or the iteration step length at any iteration time can be randomly sampled in a preset iteration step length interval, and the like. Optionally, the specified step length and the preset iteration step length interval can be set according to experience or set according to actual needs, and embodiments of the present application do not limit this.

[0134] Optionally, for the Newton method, the iteration step length is an important parameter affecting the convergence efficiency and accuracy of the algorithm, and the electronic device can also determine the iteration step length through the strong Wolfe convergence criterion, and the determined condition can be written as formula 2.4:

[0135] Formula 2.4

[0136] Wherein, the parameter t1 is usually a small value (such as t1 ∈ (0, 1)), and the parameter t2 can be greater than t1 (such as t2 ∈ (t1, 1)); optionally, t1 and t2 can be set according to experience or set according to actual needs, and embodiments of the present application do not limit this; for example, t1 can be equal to 0.001, and t2 can be equal to 0.1. Optionally, ΔX k As shown in formula 2.5:

[0137] Formula 2.5

[0138] In S205, based on the updated inversion parameter data and the current to-be-inverted parameter data, it is determined whether the convergence condition is reached.

[0139] In an implementation manner, the electronic device can use the updated inversion parameter data and the current to-be-inverted parameter data to calculate the inversion difference value, and when the inversion difference value is less than or equal to a preset convergence error, it is determined that the convergence condition is reached; when the inversion difference value is greater than the preset convergence error, it is determined that the convergence condition is not reached. Optionally, the electronic device can use formula 2.6 to calculate the inversion difference value:

[0140] Formula 2.6

[0141] Wherein, Φ can represent the inversion difference value; optionally, the preset convergence error (which can be represented as η) can be set according to experience or set according to actual needs, and embodiments of the present application do not limit this.

[0142] In another implementation, the electronic device can determine that the convergence condition is reached when the number of iterations reaches a preset iteration threshold, and / or determine that the convergence condition is not reached when the number of iterations does not reach the preset iteration threshold, and the like. The preset iteration threshold can be set according to experience or actual needs, and the present embodiment is not limited in this regard.

[0143] In S206, if the convergence condition is not reached, the updated inversion parameter data is taken as the current to-be-inverted parameter data, and the iteration of calculating the partial derivative of the target inversion function based on the current to-be-inverted parameter data is continued until the convergence condition is reached.

[0144] Based on this, the electronic device can perform k=k+1 to iteratively calculate the partial derivative of the target inversion function based on the current to-be-inverted parameter data, and the like. Optionally, in each iteration process, the electronic device can first update the target inversion function using the current to-be-inverted parameter data, that is, update the stiffness matrix using the current to-be-inverted parameter data, 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 to-be-inverted parameter data, and continuously obtain the updated inversion parameter data in the next iteration until the convergence condition is reached. It can be seen that the present embodiment can determine the current parameter value of each target parameter based on the current to-be-inverted parameter data to update the stiffness matrix, and the like.

[0145] In S207, if the convergence condition is reached, the updated inversion parameter data is taken as the target inversion result of the to-be-inverted parameter.

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

[0147] In S208, the target inversion result of each target parameter in the target parameter set is calculated based on the target inversion result of the to-be-inverted parameter, and the target parameter set includes anisotropy parameters.

[0148] Optionally, the electronic device can also determine the target stiffness matrix based on the target inversion result of each target parameter, i.e., taking the target inversion result of each target parameter as the parameter value of the corresponding target parameter to calculate the target stiffness matrix, and determine the target P-wave polarization direction coefficient and the target S-wave polarization direction coefficient based on the target stiffness matrix. Based on this, the electronic device can determine the target PP-wave diagonal parameter based on the target P-wave polarization direction coefficient, and determine the target PS-wave diagonal parameter based on the target S-wave polarization direction coefficient, and then calculate the target PP-wave reflection coefficient by using the target PP-wave diagonal parameter, the target stiffness matrix and the target P-wave polarization direction coefficient, and calculate the target PS-wave reflection coefficient by using the target PS-wave diagonal parameter, the target stiffness matrix and the target S-wave polarization direction coefficient. Optionally, the electronic device can calculate the target PP-wave reflection coefficient by using the target PP-wave diagonal parameter, the target stiffness matrix and the target P-wave polarization direction coefficient through formula 1.12, and correspondingly, calculate the target PS-wave reflection coefficient by using the target PS-wave diagonal parameter, the target stiffness matrix and the target S-wave polarization direction coefficient through formula 1.16.

[0149] In the embodiment of the present application, the target inversion result of each target parameter can be the target inversion result of each target parameter in the target region.

[0150] In the embodiment of the present application, in order to further verify the feasibility of the PP-wave and PS-wave reflection coefficient approximation and the inversion method proposed by the embodiment of the present application, on the one hand, the accuracy of the PP-wave and PS-wave reflection coefficient approximation equation is tested on a double-layer medium. For example, the upper medium can be an isotropic medium, and the lower medium can be a VTI medium with strong anisotropy characteristics. By comparing with the exact reflection coefficient (calculated by formula 1.7), it can be proved that the derived approximation equation has high accuracy, so that the calculation of the reflection coefficient can be effectively simplified while ensuring high accuracy. For example, assuming that the upper medium is an isotropic medium, the P-wave and S-wave velocities and the density are 5.5 km / s (kilometers per second), 3.5 km / s and 2.7 g / cm 3 (grams per cubic centimeter) respectively (at this time, the anisotropy parameters can all be 0), and the lower layer is a VTI medium with strong anisotropy characteristics, the P-wave and S-wave velocities and the density are 4.6 km / s, 2.5 km / s and 2.65 g / cm 3 respectively, and the anisotropy parameters are δ=0.15 and ε=0.25; in this case, the embodiment of the present application 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, as Figure 3As shown, the solid line represents the precise value of the reflection coefficient, the dotted line represents the classical weak anisotropy approximation (Rüger approximation equation), and the dashed line represents the weak contrast reflection coefficient of the PP wave proposed in this embodiment of the invention (i.e., the PP wave reflection coefficient calculated by the approximation equation proposed in this embodiment of the invention). The dashed line and the solid line basically match, proving that the PP reflection coefficient proposed in this embodiment of the invention is effective; correspondingly, as... Figure 4 As shown, the dashed line used to represent the approximate calculation of the PS wave reflection coefficient (i.e., the weak contrast reflection coefficient of the PS wave) matches the solid line better, especially at large angles, proving that the PS wave reflection coefficient proposed in this embodiment of the invention has high accuracy. Figure 3 and Figure 4 The horizontal axis represents the angle of incidence (in degrees (°)), and the vertical axis represents the reflection coefficient.

[0151] 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 well logging models 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 it with conventional anisotropic seismic prediction methods.

[0152] For example, in this embodiment of the invention, the feasibility of the inversion method was tested on a model derived from actual well logging curves of a shale reservoir. The P-wave and S-wave velocities and densities in the well logging curves can be actual well logging data, while anisotropy parameters can be calculated using a rock physics model. Optionally, seismic data can be synthesized using the accurate reflection coefficient equation (e.g., Formula 1.7), and then different levels of noise can be added to obtain seismic data with signal-to-noise ratios of 8 and 4. The joint inversion method proposed in this embodiment of the invention can then be used to predict anisotropy parameters. Based on this, Figure 5 The inversion results for seismic data with a signal-to-noise ratio of 8 are shown. The black solid line represents the true model (i.e., the true values ​​of each target parameter in the corresponding test area, such as those calculated using a rock physics model), the gray solid line represents the initial model (e.g., obtained based on experience or testing), and the dashed line represents the target inversion results. Figure 5 As can be seen, the target inversion results match the actual model well, meaning the accuracy of the target inversion results for each target parameter (such as the target inversion results for P-wave velocity) is high; furthermore, Figure 6 The inversion results for seismic data with a signal-to-noise ratio (SNR) of 4 show that the inversion results for each target parameter are close to the true values, indicating high accuracy. It should be understood that compared to the inversion results with an SNR of 8, the inversion accuracy of the results with an SNR of 4 is slightly lower, but it still closely matches the true model, especially in areas with high anisotropy, demonstrating the feasibility and effectiveness of the proposed method. Figure 5 and Figure 6The ordinate is time and the abscissa is amplitude.

[0153] In another aspect, the joint inversion method proposed by the embodiments of the present application is also applied to a practical shale work area, and by comparison with the inversion method based on the weak anisotropy approximation, the reliability of the joint inversion method proposed by the embodiments of the present application can be effectively verified. Among them, Figure 7 The anisotropy parameter δ inverted by the conventional inversion method (such as the Rüger inversion method) is shown by the black line, and the anisotropy parameter δ calculated by the rock physics model is shown by the black line (the black line is used to indicate that the parameter value indicated by the curve located on the right (i.e. the right of the dotted line) is higher, and the parameter value indicated by the curve located on the left is lower, that is, the higher the black line deviates to the right, the larger the parameter value indicated by the black line), it can be seen that there is a certain difference between the result obtained by the conventional method and the logging curve (i.e. the black line), wherein the highest value on the logging curve indicates the location of the shale reservoir, and the color of the conventional inversion method at the highest value on the logging curve 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 8 The anisotropy parameter δ obtained by the joint inversion method proposed by the embodiments of the present application is shown by the black line, and the inversion result is in good agreement with the logging curve, and the location of the shale reservoir is clearly indicated, that is, the target inversion result of the anisotropy parameter δ reflects that the brightness at the highest value on the black line is the brightest (i.e. 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 that 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, Figure 9 and Figure 10 The anisotropy parameters ε inverted by the conventional method and the joint inversion method proposed by the embodiments of the present application are shown by the black line and the black line, respectively; compared with the conventional method, the target inversion result of the anisotropy parameter ε obtained by the joint inversion method proposed by the embodiments of the present application is in better agreement with the logging curve, which proves that the inversion proposed by the embodiments of the present application can effectively improve the stability and precision of the anisotropy parameter prediction. Among them, Figures 7-10 The ordinate can be time and the abscissa can be seismic trace number.

[0154] The embodiment of the present application can construct a forward operator based on the initial parameter data after obtaining the initial parameter data of the target region and obtaining the PP wave seismic data and the PS wave seismic data of the target region, 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 to-be-inverted parameter data of the to-be-inverted parameter can be determined, and the partial derivative of the target inversion function can be calculated based on the current to-be-inverted parameter data; and the to-be-inverted parameter can be updated based on the partial derivative of the target inversion function to obtain the updated to-be-inverted parameter data of the to-be-inverted parameter. Based on this, whether the convergence condition is reached can be judged based on the updated to-be-inverted parameter data and the current to-be-inverted parameter data; if the convergence condition is not reached, the updated to-be-inverted parameter data is taken as the current to-be-inverted parameter data, and the iteration is continued to perform the calculation of the partial derivative of the target inversion function based on the current to-be-inverted parameter data until the convergence condition is reached; if the convergence condition is reached, the updated to-be-inverted parameter data is taken as the target inversion result of the to-be-inverted parameter. Further, the target inversion result of each target parameter in the target parameter set can be calculated based on the target inversion result of the to-be-inverted parameter, and the target parameter set includes anisotropy parameters. It can be seen that the embodiment of the present application can realize joint inversion through the PP wave seismic data and the PS wave seismic data, and can effectively improve the stability and precision of anisotropy parameter prediction.

[0155] Based on the description of the related embodiments of the above joint inversion method, the embodiment of the present application further proposes a joint inversion device, which can be a computer program (including program code) running in an electronic device; as shown in the figure, the joint inversion device can include an acquisition unit 1101 and a processing unit 1102. The joint inversion device can perform the joint inversion method as shown in the figure, that is, the joint inversion device can run the above units: Figure 11 Figure 1 or Figure 2 The joint inversion device can perform the joint inversion method as shown in the figure, that is, the joint inversion device can run the above units:

[0156] The acquisition unit 1101 is used to acquire the initial parameter data of the target region, and acquire the PP wave seismic data and the PS wave seismic data of the target region, and the initial parameter data includes upper medium parameter data and lower medium parameter data;

[0157] The processing unit 1102 is used to 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;

[0158] The processing unit 1102 is further used to determine the target inversion result of the to-be-inverted parameter based on the target inversion function;

[0159] ​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, the target parameter set including an anisotropy parameter.

[0160] In an embodiment, the forward operator includes a first forward operator and a second forward operator; and the processing unit 1102, when constructing the forward operator based on the initial parameter data, can be specifically configured to:

[0161] determine a stiffness matrix based on the initial parameter data, the stiffness matrix including a plurality of stiffness elements, an initial element value of a stiffness element being 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 including the initial value of each target parameter in the target parameter set, the initial value of a target parameter including an upper initial value of the corresponding target parameter in the upper medium parameter data and / or a lower initial value of the corresponding target parameter in the lower medium parameter data;

[0162] calculate a plurality of PP wave parameter perturbation coefficients based on the stiffness matrix; and calculate a plurality of PS wave parameter perturbation coefficients based on the stiffness matrix; wherein the plurality of PP wave parameter perturbation coefficients include 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 PS wave parameter perturbation coefficients corresponding to each stiffness element;

[0163] construct the first forward operator based on the plurality of PP wave parameter perturbation coefficients; and construct the second forward operator based on the plurality of PS wave parameter perturbation coefficients, so as to construct the forward operator.

[0164] In another embodiment, the processing unit 1102, when calculating the plurality of PP wave parameter perturbation coefficients based on the stiffness matrix, can be specifically configured to:

[0165] determine a plurality of intermediate coefficients based on the stiffness matrix; and calculate a PP wave diagonal parameter based on the plurality of intermediate coefficients, a P-wave polarization direction coefficient, and a P-wave vertical slowness;

[0166] calculate the plurality of PP wave parameter perturbation coefficients based on the PP wave diagonal parameter, the stiffness matrix, and the plurality of intermediate coefficients;

[0167] The processing unit 1102, when calculating the plurality of PS wave parameter perturbation coefficients based on the stiffness matrix, can be specifically configured to:

[0168] calculate a PS wave diagonal parameter based on the plurality of intermediate coefficients, a S-wave polarization direction coefficient, and a S-wave vertical slowness;

[0169] Based on the PS wave diagonal parameter, the stiffness matrix and the plurality of intermediate coefficients, a plurality of PS wave parameter perturbation coefficients are calculated.

[0170] In another implementation, the forward operator includes a first forward operator and a second forward operator; and the processing unit 1102, when constructing the target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data, can be specifically configured to:

[0171] determine a stiffness matrix and an anisotropy parameter relationship constraint term;

[0172] construct a PP wave seismic data error term based on the first forward operator and the PP wave seismic data, and construct a PS wave seismic data error term based on the second forward operator and the PS wave seismic data;

[0173] construct the target inversion function by using the PP wave seismic data error term, the PS wave seismic data error term, and the stiffness matrix and the anisotropy parameter relationship constraint term.

[0174] In another implementation, the processing unit 1102, when determining the target inversion result of the to-be-inverted parameter based on the target inversion function, can be specifically configured to:

[0175] determine current to-be-inverted parameter data of the to-be-inverted parameter, and calculate a partial derivative of the target inversion function based on the current to-be-inverted parameter data;

[0176] update the to-be-inverted parameter based on the partial derivative of the target inversion function to obtain updated inversion parameter data of the to-be-inverted parameter;

[0177] determine whether a convergence condition is reached based on the updated inversion parameter data and the current to-be-inverted parameter data;

[0178] if the convergence condition is not reached, the updated inversion parameter data is taken as the current to-be-inverted parameter data, and the iteration of calculating the partial derivative of the target inversion function based on the current to-be-inverted parameter data is continued until the convergence condition is reached;

[0179] if the convergence condition is reached, the updated inversion parameter data is taken as the target inversion result of the to-be-inverted parameter.

[0180] In another implementation, the processing unit 1102, when calculating the target inversion result of each target parameter in a target parameter set based on the target inversion result of the to-be-inverted parameter, can be specifically configured to:

[0181] determine a calculation mode corresponding to each target parameter in the target parameter set based on a nonlinear relationship between a stiffness matrix and the target parameter set, respectively.

[0182] The target inversion result of each target parameter is calculated according to the calculation mode corresponding to the target parameter and the target inversion result of the parameter to be inverted.

[0183] In another implementation, the processing unit 1102 can also be configured to:

[0184] The target stiffness matrix is determined based on the target inversion result of each target parameter;

[0185] The target P-wave polarization direction coefficient and the target S-wave polarization direction coefficient are determined based on the target stiffness matrix;

[0186] The target PP-wave diagonal parameter is determined based on the target P-wave polarization direction coefficient, and the target PS-wave diagonal parameter is determined based on the target S-wave polarization direction coefficient;

[0187] The target PP-wave reflection coefficient is calculated using the target PP-wave diagonal parameter, the target stiffness matrix, and the target P-wave polarization direction coefficient, and the target PS-wave reflection coefficient is calculated using the target PS-wave diagonal parameter, the target stiffness matrix, and the target S-wave polarization direction coefficient.

[0188] According to one embodiment of the present application, Figure 11 The units in the joint inversion device shown can be combined into one or several other units to constitute, respectively or all, or some of the units can be further split into a plurality of units with smaller functions to constitute, which can achieve the same operation without affecting the implementation of the technical effects of the embodiments of the present application. The above units are divided based on logical functions, and in actual application, the functions 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 application, any joint inversion device can also include other units, and in actual application, these functions can also be assisted by other units, and can be implemented by multiple units.

[0189] According to another embodiment of the present application, the joint inversion device as shown in Figure 1 or Figure 2 the corresponding method shown in can be constructed, and the joint inversion method of the embodiments of the present application can be implemented. The computer program (including program code) involved in each step of the corresponding method shown in Figure 11 the joint inversion device shown in can be recorded on, for example, a computer storage medium, and loaded into the above-mentioned electronic device through the computer storage medium, and run therein.

[0190] Based on the description of the method embodiments and the device embodiments, the exemplary embodiments of the present application further provide an electronic device, comprising: at least one processor; and a memory connected with the at least one processor in communication. The memory stores a computer program capable of being executed by the at least one processor, and the computer program, when executed by the at least one processor, is configured to cause the electronic device to perform the method according to the embodiments of the present application.

[0191] The exemplary embodiments of the present application further 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 configured to cause the computer to perform the method according to the embodiments of the present application.

[0192] The exemplary embodiments of the present application further provide a computer program product comprising a computer program, wherein the computer program, when executed by a processor of a computer, is configured to cause the computer to perform the method according to the embodiments of the present application.

[0193] Reference Figure 12 will now be described, which is an example of a hardware device that can be applied to various aspects of the present application. The electronic device is intended to represent various forms of digital electronic computing devices such as laptops, desktops, workstations, personal digital assistants, servers, blade servers, mainframes, and other appropriate computing devices. The electronic device can also represent various forms of mobile devices such as personal digital processing, cellular telephones, smart phones, wearable devices, and other similar computing devices. The components, their connections, and their functions, as shown in the figures, and described in this document, are meant only to be examples and are not intended to limit the implementations of the application described and / or claimed in this document.

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

[0195] The various 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 capable of inputting information to the electronic device 1200, and can receive inputted digital or character information, as well as generate key signal inputs related to user settings and / or function controls of the electronic device. The output unit 1207 can be any type of device capable of presenting 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 magnetic disk, an optical disk. The communication unit 1209 allows the electronic device 1200 to exchange information / data with other devices over 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.

[0196] The computing unit 1201 can be various general and / or special purpose processing components having 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 specialized artificial intelligence (AI) computing chips, various computing units running machine learning model algorithms, a digital signal processor (DSP), and any appropriate processor, controller, microcontroller, etc. The computing unit 1201 performs various methods and processes described above. For example, in some embodiments, the joint inversion method can be implemented as a computer software program tangibly embodied in a machine-readable medium, such as the storage unit 1208. In some embodiments, part or all of the computer program can be loaded and / or installed onto the electronic device 1200 via the ROM 1202 and / or the communication unit 1209. In some embodiments, the computing unit 1201 can be configured to perform the joint inversion method by any other appropriate means, such as by means of firmware.

[0197] Program code for carrying out the methods of the present application can be written in any combination of one or more programming languages. The program code can be provided to a processor or controller of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the program code, when executed by the processor or controller, causes the machine to perform the functions / operations specified in the flow diagrams and / or block diagrams. The program code can execute entirely on a machine, partly on a machine, as a stand-alone software package, partly on a machine and partly on a remote machine or entirely on a remote machine or server.

[0198] In the context of the present application, a machine-readable medium can be a tangible medium that contains or stores a program for use by or in connection with an instruction execution system, apparatus, or device. The machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. A machine-readable medium can include but is not limited to an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples of the machine-readable storage medium would include an electrical connection based on one or more wires, 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 disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.

[0199] As used in the present application, the terms "machine-readable medium" and "computer-readable medium" refer to any computer program product, apparatus and / or device (e.g., magnetic discs, optical disks, memory, Programmable Logic Devices (PLDs)) used to provide 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 that can be used to provide machine instructions and / or data to a programmable processor.

[0200] To provide for interaction with a user, the systems and techniques described here 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 a pointing device (e.g., a mouse or a trackball) by which the user can provide input to the computer. Other kinds of devices can be used to provide for interaction with a user as well; for example, feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user can be received in any form, including acoustic, speech, or tactile input.

[0201] Also, it is to be understood that the above-revealed embodiments are merely exemplary of the application and should not be used to limit the scope of the present application, which is defined by the claims. Any equivalent variation of the above-described arrangements falls within the scope of the present application.

Claims

1. A joint inversion method, characterized in that, include: Acquire initial parameter data for the target area, and acquire PP wave seismic data and PS wave seismic data for the target area. The initial parameter data includes upper medium parameter data and lower medium parameter data. Based on the initial parameter data, a forward modeling operator is constructed, including: determining a stiffness matrix based on the initial parameter data, the stiffness matrix including multiple stiffness elements, the initial element value of a stiffness element being determined based on the initial value of each target parameter in at least one target parameter in the target parameter set, the initial parameter data including the initial values ​​of each target parameter in the target parameter set, the initial value of a target parameter including the upper-level initial value of the corresponding target parameter in the upper-level medium parameter data and / or the lower-level initial value in the lower-level medium parameter data; wherein, the initial element value of a stiffness element being determined based on the initial value of each target parameter in at least one target parameter in the target parameter set means that: the initial element value of a stiffness element is determined based on the initial parameter value of each target parameter in at least one target parameter in the target parameter set, and the initial parameter value of a target parameter is the average between the upper-level initial value and the lower-level initial value of the corresponding target parameter; Based on the stiffness matrix, multiple PP wave parameter perturbation coefficients are calculated; and based on the stiffness matrix, multiple PS wave parameter perturbation coefficients are calculated; wherein, the multiple PP wave parameter perturbation coefficients include the PP wave parameter perturbation coefficients corresponding to each stiffness element in the multiple stiffness elements, and the multiple PS wave parameter perturbation coefficients include the PS wave parameter perturbation coefficients corresponding to each stiffness element. Based on the perturbation coefficients of the multiple PP wave parameters, a first forward modeling operator is constructed; and based on the perturbation coefficients of the multiple PS wave parameters, a second forward modeling operator is constructed to realize the construction of the forward modeling operator; and based on the forward modeling operator, the PP wave seismic data, and the PS wave seismic data, a target inversion function is constructed, including: determining the stiffness matrix and anisotropic parameter relationship constraint terms; Based on the perturbation coefficients of the multiple PP wave parameters, a first forward modeling operator is constructed; and based on the perturbation coefficients of the multiple PS wave parameters, a second forward modeling operator is constructed to realize the construction of the forward modeling operator; and based on the forward modeling operator, the PP wave seismic data, and the PS wave seismic data, a target inversion function is constructed, including: determining the stiffness matrix and anisotropic parameter relationship constraint terms; Based on the first forward modeling operator and the PP wave seismic data, a PP wave seismic data error term is constructed; and based on the second forward modeling operator and the PS wave seismic data, a PS wave seismic data error term is constructed. The target inversion function is constructed using the error terms of the PP wave seismic data, the error terms of the PS wave seismic data, and the constraint terms on the relationship between the stiffness matrix and the anisotropic parameters. Based on the target inversion function, the target inversion result of the parameters to be inverted is determined; based on the target inversion result of the parameters to be inverted, the target inversion result of each target parameter in the target parameter set is calculated, wherein the target parameter set includes anisotropic parameters; wherein, the target inversion function is: ； Where, d p The PP wave seismic data is given, G1 is the first forward modeling operator, X is the parameter to be inverted, and d is the parameter to be inverted. s For the PS wave seismic data, G2 is the second forward modeling operator, α is the constraint weight, and the function g represents the function that calculates the anisotropy parameters using the stiffness matrix; The step of constructing the first forward modeling operator based on the plurality of PP wave parameter perturbation coefficients includes: extracting a first wavelet from the PP wave seismic data, including: superimposing seismic data from each angle trace in the PP wave seismic data to obtain PP wave post-stack seismic data, so as to extract the first wavelet from the PP wave post-stack seismic data; and constructing the first forward modeling operator based on the first wavelet and the plurality of PP wave parameter perturbation coefficients; the step of constructing the second forward modeling operator based on the plurality of PS wave parameter perturbation coefficients includes: extracting a second wavelet from the PS wave seismic data, including: superimposing seismic data from each angle trace in the PS wave seismic data to obtain PS wave post-stack seismic data, so as to extract the second wavelet from the PS wave post-stack seismic data; and constructing the second forward modeling operator based on the second wavelet and the plurality of PS wave parameter perturbation coefficients; The step of determining the target inversion result of the parameters to be inverted based on the target inversion function includes: determining the current data of the parameters to be inverted, and calculating the partial derivative of the target inversion function based on the current data of the parameters to be inverted; Based on the partial derivative of the target inversion function, the parameters to be inverted are updated to obtain the updated inversion parameter data of the parameters to be inverted. Based on the updated inversion parameter data and the current parameter data to be inverted, it is determined whether the convergence condition has been met; if the convergence condition has not been met, the updated inversion parameter data is used as the current parameter data to be inverted, and the process of calculating the partial derivative of the target inversion function based on the current parameter data to be inverted continues iteratively 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 parameters to be inverted; wherein, the partial derivatives of the target inversion function include the first-order partial derivatives and the second-order partial derivatives of the target inversion function, and the formulas for calculating the first-order partial derivatives and the second-order partial derivatives are as follows: ； Among them, ▽f and ▽ 2 f represents the first-order partial derivative and the second-order partial derivative, respectively; J and H are the Jacobian matrix and Hessian matrix of the function g, respectively; G P G represents the longitudinal wave polarization direction coefficient. S Represents the transverse wave polarization direction coefficient; calculates the perturbation coefficients of multiple PP wave parameters, including: Based on the stiffness matrix, multiple intermediate coefficients are determined; and based on the multiple intermediate coefficients, the longitudinal wave polarization direction coefficient, and the longitudinal wave vertical slowness, the PP wave diagonal parameters are calculated. Based on the PP wave diagonal parameters, the stiffness matrix, and the multiple intermediate coefficients, calculate multiple PP wave parameter perturbation coefficients; The calculation of multiple PS wave parameter perturbation coefficients based on the stiffness matrix includes: Based on the aforementioned intermediate coefficients, transverse wave polarization direction coefficients, and transverse wave vertical slowness, the PS wave diagonal parameters are calculated. Based on the PS wave diagonal parameters, the stiffness matrix, and the multiple intermediate coefficients, calculate multiple PS wave parameter perturbation coefficients; The calculation of the target inversion result for each target parameter in the target parameter set based on the target inversion result of the parameters to be inverted includes: Based on the nonlinear relationship between the stiffness matrix and the target parameter set, the calculation method for each target parameter in the target parameter set is determined respectively; 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; Based on the target inversion results of the aforementioned target parameters, the target stiffness matrix is ​​determined. Based on the target stiffness matrix, the target longitudinal wave polarization direction coefficient and the target transverse wave polarization direction coefficient are determined; The target PP wave diagonal parameters are determined based on the target longitudinal wave polarization direction coefficient, and the target PS wave diagonal parameters are determined based on the target transverse wave polarization direction coefficient. The target PP wave reflection coefficient is calculated using the target PP wave diagonal parameters, the target stiffness matrix, and the target longitudinal wave polarization direction coefficient; and the target PS wave reflection coefficient is calculated using the target PS wave diagonal parameters, the target stiffness matrix, and the target transverse wave polarization direction coefficient.

Images