Stepwise Inversion Method for Orthotropic Media Based on the Second Derivative of Elastic Impedance

Abstract

The present invention relates to the field of oil and gas exploration and development, and particularly to a step-by-step inversion method for orthotropic media based on the second derivative of elastic impedance. The method includes the following steps: First, estimate the corresponding azimuthal elastic impedance data from the azimuthal seismic data of partial angle stacking; Second, calculate the ratio of the elastic impedance data of two different azimuths, and solve the normal weakness parameter and tangential weakness parameter of the fracture; Third, based on the azimuthal elastic impedance equation of the orthotropic medium with parameter combination and mathematical approximation, calculate the elastic impedance data containing only the incident angle information, invert attributes A, B, C, and D, and obtain the P-wave anisotropy parameter by using the ratio of attributes A and C. The inversion method proposed by the present invention improves the prediction accuracy of model parameters by introducing the second derivative of elastic impedance with respect to model parameters, and reduces the number of model parameters to be inverted by using the step-by-step inversion strategy, thereby improving the reliability of model parameter prediction.

Description

Technical Field

[0001] The present invention relates to the field of oil and gas exploration and development, and particularly to a step-by-step inversion method for orthotropic media based on the second derivative of elastic impedance. Background Art

[0002] Fractures control the permeability of underground reservoirs and fluid migration, and strongly affect rock properties. Seismic inversion is an important tool for obtaining underground fracture information and realizing reservoir spatial description. Therefore, using seismic inversion methods for fracture detection plays a crucial role in well location optimization, reservoir exploration management, carbon dioxide storage and other fields.

[0003] In recent years, a large number of field outcrops and logging data have confirmed that rocks with vertical and near-vertical fractures develop in a vertically transverse isotropic (VTI) background (such as periodic thin interbeds, horizontal fractures, shale, etc.) that is prevalent underground. Such rocks can usually be equivalently regarded as orthotropic (ORT) media. The azimuthal seismic reflection coefficient equation is a bridge connecting the anisotropic characteristics of reservoirs and the macroscopic seismic response, laying a solid mathematical and physical foundation for the construction of seismic forward and inverse operators.

[0004] However, the existing PP-wave azimuthal seismic reflection coefficient equation for describing orthotropic media contains many model parameters, and the weight coefficients between the model parameters vary greatly, making it difficult to stably solve these parameters simultaneously during seismic inversion, thereby affecting the exploration and development and fine description requirements of orthotropic media reservoirs.

[0005] Specifically, under the assumption of long seismic wavelengths, rocks with a set of vertical fractures in a VTI background can be equivalent to ORT media. Under the assumption of weak anisotropy, elastic parameters are usually used in the industry to describe isotropic backgrounds, and anisotropy parameters and fracture parameters are used to describe the strength of anisotropy in the VTI background and the degree of development of vertical fractures, respectively. At this time, the reflection coefficient equation describing the ORT medium usually contains 8 unknown parameters. Even if the mutual coupling effect between vertical fractures is ignored, the equation still contains 7 model parameters, and the weight coefficients between the model parameters vary greatly, making it difficult to stably solve these parameters simultaneously during seismic inversion. Existing technologies mostly slow down the uncertainty of equation solving based on parameter reduction and multi-parameter inversion strategies. (1) In terms of the parameter reduction strategy: The existing multi-parameter reduction strategies for reflection coefficient equations usually discard large incident angle information during the parameter combination process, resulting in limited accuracy in large incident angle seismic inversion applications. In addition, it is difficult for existing technologies to directly predict anisotropy parameters related to the VTI background from the combined attribute parameters. (2) In terms of the multi-parameter inversion strategy: The existing step-by-step inversion methods can effectively reduce the number of model parameters and reduce the uncertainty of multi-parameter inversion. However, the variation of amplitude with incident angle and azimuth (Amplitude variation angle and azimuth, abbreviated as AVAZ) is sensitive to noise, resulting in still controversial applicability and reliability of AVAZ inversion in the presence of noise. Although the existing azimuthal elastic impedance inversion methods can improve the noise resistance of seismic inversion, when using the elastic impedance inversion model parameters, the elastic impedance is usually logarithmized, and then a nonlinear relationship between the logarithmic elastic impedance and the logarithmic elastic parameters is established. Its theoretical basis is the first-order linear approximation of the reflection coefficient, and the inversion accuracy only depends on the weight coefficients of the parameters in the first-order linear approximation formula to a certain extent. In summary, existing technologies are difficult to meet the actual exploration, development and fine description requirements of orthogonally fractured reservoirs. Summary of the Invention

[0006] The technical problem to be solved by the present invention is to overcome the defects existing in the prior art, and propose an azimuthal elastic impedance equation for orthotropic media based on parameter combination and mathematical approximation. When inverting, the prediction accuracy of model parameters is improved by introducing the second derivative of elastic impedance with respect to model parameters, and the step-by-step inversion strategy is used to reduce the number of model parameters to be inverted, thereby improving the reliability of model parameter prediction.

[0007] The technical solution adopted by the present invention to solve the above technical problem is: A step-by-step inversion method for orthotropic media based on the second derivative of elastic impedance, which includes the following steps:

[0008] In the first step, estimate the corresponding azimuthal elastic impedance data from the azimuthal seismic data stacked at partial angles;

[0009] In the second step, calculate the ratio of the elastic impedance data in two different orientations, and solve for the normal weakness parameter and tangential weakness parameter of the fracture;

[0010] In the third step, based on the azimuthal elastic impedance equation of an orthotropic medium with parameter combinations and mathematical approximations, calculate the elastic impedance data that only contains incident angle information, invert attributes A, B, C, and D, and calculate the P-wave anisotropy parameter ε using the ratio of attributes A and C;

[0011] The azimuthal elastic impedance equation of an orthotropic medium based on parameter combinations and mathematical approximations is as follows:

[0012]

[0013] In Equation (5), the symbol "exp" represents the exponential operation, and the subscript symbol "0" represents the constant elastic parameters of the background medium.

[0014]

[0015] In Equations (5) and (6), AEI represents the azimuthal elastic impedance, θ represents the incident angle, represents the azimuth angle, M and μ represent the P-wave modulus and S-wave modulus of the isotropic background, ρ represents the density, g = μ / M; Attribute A represents the P-wave modulus; Attribute B represents the anisotropic shear modulus, and when the anisotropy parameter is 0, Attribute B represents the isotropic background shear modulus; Attribute C represents the anisotropic P-wave modulus; Attribute D represents the density; δ N represents the normal weakness parameter of the fracture; δ T represents the tangential weakness parameter of the fracture; ε represents the P-wave anisotropy parameter.

[0016] Preferably, in the first step:

[0017] For seismic data containing i sampling points, the relationship between the seismic data and the azimuthal elastic impedance data is represented in matrix form as:

[0018] S PP = G1X (8),

[0019] In Equation (8),

[0020]

[0021] In Equation (9), j and k represent the number of incident angles and azimuth angles respectively, s i represents the seismic record of the i-th sampling point, lnAEI represents the logarithm of the azimuthal elastic impedance, W represents the wavelet matrix, D o represents the difference operation matrix;

[0022] When the initial model X of azimuthal elastic impedance is given pri After that, Equation (8) is solved by the least squares algorithm based on model parameters and is expressed as:

[0023]

[0024] In Equation (10), represents the inverse of the covariance matrix of model parameters, σ represents the damping factor related to the signal-to-noise ratio, and Equation (10) is used for solution to estimate azimuthal elastic impedance data from azimuthal seismic data with partial angle stacking.

[0025] Preferably, in the second step:

[0026] For two different observation azimuths and azimuth The ratio of elastic impedance is expressed as:

[0027]

[0028] In Equation (10), DEI represents the ratio of elastic impedance data at two different azimuths,

[0029]

[0030] Equation (11) is expressed in matrix form as:

[0031] d f = G2m f (13),

[0032] In Equation (13), G2 represents the non-linear forward operator related to the incident angle and azimuth angle, and

[0033]

[0034] The model parameters of Equation (13) are solved by the Newton method and expressed as:

[0035] m l+1 = m l + β l Δm l (15),

[0036] In Equation (15), β l represents the step size of each iteration, which is obtained by the line search method, m l+1 represents the predicted result of model parameters after the l-th iteration, m l represents the initial value of model parameters at the l-th iteration, Δm l respectively represent the perturbation amount of model parameters at the l-th iteration and are expressed as:

[0037] Δm = -H -1Y(16),

[0038] In Equation (16), Y and H respectively represent the first-order derivative and second-order derivative of the model observation data vector d mod (m) with respect to the model parameter m, and the approximate solution is:

[0039]

[0040] [[ID=eleven]]In Equation (16), d input (m) represents the input actual observation data vector, and the symbol "diag" represents a diagonal matrix;

[0041] According to Equations (11) and (13), in Equation (17) is expressed as:

[0042]

[0043] In Equation (18),

[0044]

[0045] Substitute Equation (18) into Equations (15)-(17) to obtain the normal weakness parameter and tangential weakness parameter of the crack.

[0046] Preferably, in the third step:

[0047] The elastic impedance equation related only to the incident angle is expressed as:

[0048] )]]

[0049] Equation (20) is expressed in matrix form as:

[0050] d v = G3m v (21),

[0051] In Equation (20), G3 represents a non-linear forward operator related only to the incident angle, and

[0052]

[0053] According to Equations (20) and (21), in Equation (17) is expressed as:

[0054]

[0055] In Equation (23),

[0056]

[0057] Substitute Equation (23) into Equations (15)-(17) to obtain the parameters of Property A, Property B, Property C, and Property D;

[0058] Use Equation to calculate the P-wave anisotropy parameter ε.

[0059] The beneficial technical effects brought by this invention:

[0060] The azimuthal elastic impedance equation of orthotropic media based on parameter combination and mathematical approximation proposed by this invention contains fewer model parameters, each of which has a strict physical meaning, and the P-wave anisotropy parameter related to the VTI background can be directly calculated through the ratio of Property A and Property C.

[0061] Directly calculate the second derivative of the model parameters according to the azimuthal elastic impedance equation proposed by this invention. The calculation process is simple, and it can improve the prediction accuracy of the model parameters.

[0062] The method of this invention reduces the number of model parameters in each inversion by introducing a step-by-step inversion strategy, reduces the uncertainty of the inversion of orthotropic media, and thus improves the reliability of the prediction of the model parameters in the equation. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is a flowchart of an embodiment of this invention.

[0064] Figure 2 is a schematic diagram of a synthetic azimuthal seismic gather. Among them, (a) is a schematic diagram of a synthetic azimuthal seismic gather with a signal-to-noise ratio of 10:1, and (b) is a schematic diagram of a synthetic azimuthal seismic gather with a signal-to-noise ratio of 4:1.

[0065] Figure 3 is a schematic diagram of the azimuthal elastic impedance inversion result with a signal-to-noise ratio of 10:1. Among them, (a) is a schematic diagram of the elastic impedance inversion result at different incident angles with an azimuth angle of 0°, and (b) is a schematic diagram of the elastic impedance inversion result at different incident angles with an azimuth angle of 90°.

[0066] Figure 4 is a schematic diagram of the azimuthal elastic impedance inversion result with a signal-to-noise ratio of 4:1. Among them, (a) is a schematic diagram of the elastic impedance inversion result at different incident angles with an azimuth angle of 0°, and (b) is a schematic diagram of the elastic impedance inversion result at different incident angles with an azimuth angle of 90°.

[0067] Figure 5Schematic diagram of the inversion results of the model parameters at a signal-to-noise ratio of 10:1. Among them, (a) is the schematic diagram of the inversion result of property A, (b) is the schematic diagram of the inversion result of property B, (c) is the schematic diagram of the inversion result of property C, (d) is the schematic diagram of the inversion result of property D, (e) is the schematic diagram of the inversion result of the crack normal weakness parameter, (f) is the schematic diagram of the inversion result of the crack tangential weakness parameter, and (g) is the schematic diagram of the inversion result of the P-wave anisotropy parameter.

[0068] Figure 6 Schematic diagram of the inversion results of the model parameters at a signal-to-noise ratio of 4:1. Among them, (a) is the schematic diagram of the inversion result of property A, (b) is the schematic diagram of the inversion result of property B, (c) is the schematic diagram of the inversion result of property C, (d) is the schematic diagram of the inversion result of property D, (e) is the schematic diagram of the inversion result of the crack normal weakness parameter, (f) is the schematic diagram of the inversion result of the crack tangential weakness parameter, and (g) is the schematic diagram of the inversion result of the P-wave anisotropy parameter.

[0069] Figure 7 Schematic diagram of the inversion result profile. Among them, (a) is the schematic diagram of the inversion result profile of property A, (b) is the schematic diagram of the inversion result profile of property B, (c) is the schematic diagram of the inversion result profile of property C, (d) is the schematic diagram of the inversion result profile of property D, (e) is the schematic diagram of the inversion result profile of the crack normal weakness parameter, (f) is the schematic diagram of the inversion result profile of the crack tangential weakness parameter, and (g) is the schematic diagram of the inversion result profile of the P-wave anisotropy parameter. Detailed implementation method

[0070] The detailed description and technical content of the present invention are described below in conjunction with the accompanying drawings. However, the accompanying drawings are only provided for reference and explanation, and are not used to limit the present invention.

[0071] In the present invention, the derivation of the azimuthal elastic impedance equation for the orthotropic ORT medium is mainly to construct the azimuthal elastic impedance equation for the PP wave of the orthotropic ORT medium through parameter combination and mathematical approximation methods.

[0072] Under the weak anisotropy approximation and the assumption of small perturbations of elastic parameters on both sides of the reflection interface, the azimuthal seismic reflection coefficient equation for an orthotropic medium composed of rocks with a set of vertical cracks developed in the VTI background can be expressed as:

[0073]

[0074] Among them, R PP represents the azimuthal seismic reflection coefficient; θ represents the incident angle, Denote the azimuth angle; M and μ denote the longitudinal wave modulus and the transverse wave modulus of the isotropic background, ρ denotes the density; ε and δ denote the anisotropic parameters of the VTI background, where ε is used to measure the longitudinal wave anisotropy and is also called the longitudinal wave anisotropic parameter; δ N , δ V and δ H respectively denote the normal weakness parameter, the vertical tangential weakness parameter and the horizontal tangential weakness parameter of the fracture; the symbol “—” at the top of the model parameter denotes the average value of the model parameters of the upper and lower layer media; the symbol “Δ” denotes the difference of the model parameters of the upper and lower layer media; g = μ / M.

[0075] Assume that the fracture is a rotation-invariant fracture (that is, the fractures are arranged in parallel alignment and there is no cross-coupling effect between the vertical fractures), then there is δ V = δ H = δ T , where δ T denotes the tangential weakness parameter of the fracture. In addition, note that, sin 2 θ tan 2 θ = tan 2 θ - sin 2 θ, so equation (1) can be re-expressed as:

[0076]

[0077] where,

[0078]

[0079] The approximate relationship between the azimuthal seismic reflection coefficient and the azimuthal elastic impedance can be expressed as:

[0080] <U+

[0081] where, AEI denotes the azimuthal elastic impedance.

[0082] Furthermore, considering the mathematical approximation Δx / x ≈ Δlnx (basic operation in the industry), and combining equations (2), (3) and (4) for mathematical operations, and taking the integral of the operation results, and then combining and simplifying the model parameters, the normalized azimuthal elastic impedance equation of the orthotropic ORT medium can be obtained:

[0083]

[0084] Equation (5) is the azimuthal elastic impedance equation of the orthotropic medium based on parameter combination and mathematical approximation of the present invention, where the symbol “exp” represents the exponential operation, and the subscript symbol “0” represents the elastic parameters of the constant background medium, which can usually be obtained by averaging the corresponding logging curves (Whitcombe, 2002), and:

[0085]

[0086] Equation (5) only contains six parameters, and each property parameter has a clear physical meaning: Property A represents the longitudinal wave modulus; Property B represents the anisotropic shear modulus. When the anisotropy parameter is 0, Property B can be expressed as the isotropic background shear modulus; Property C represents the anisotropic longitudinal wave modulus; Property D represents the density; δ N and δ T represent the normal weakness parameter and the tangential weakness parameter of the fracture, which can be used to estimate the fracture density, calculate the fluid factor, etc. After obtaining Property A and Property C, the longitudinal wave anisotropy parameter ε directly related to the VTI background can also be obtained by their ratio, that is:

[0087]

[0088] As Figure 1 shown, a step-by-step inversion method for orthotropic media based on the second derivative of elastic impedance according to an embodiment of the present invention includes the following steps:

[0089] First step, estimate the corresponding azimuth elastic impedance data (AEI) from the azimuth seismic data of partial angle stacking;

[0090] Second step, calculate the ratio of the elastic impedance data of two different azimuths, that is, the elastic impedance data (DEI) that only contains azimuth difference information, and solve the normal weakness parameter and the tangential weakness parameter of the fracture;

[0091] Third step, based on the azimuth elastic impedance equation of the orthotropic medium with parameter combination and mathematical approximation, calculate the elastic impedance data (EI) that only contains incident angle information, invert Property A, Property B, Property C, Property D, and use the ratio of Property A and Property C to obtain the longitudinal wave anisotropy parameter ε.

[0092] First step: Estimate the azimuth elastic impedance data from the azimuth seismic data of partial angle stacking.

[0093] For seismic data containing i sampling points, the relationship between the seismic data and the azimuth elastic impedance data can be expressed in matrix form as:

[0094] S PP = G1X(8)

[0095] Where

[0096]

[0097] In the formula, j and k respectively represent the number of incident angles and azimuth angles. s iDenote the seismic record of the \(i\)-th sampling point, where the superscript "T" represents the transpose of a matrix. \(\ln AEI\) represents the logarithm of the azimuthal elastic impedance \(AEI\). \(W\) represents the wavelet matrix, and \(D\) o represents the difference operation matrix.

[0098] When the initial model \(X\) of the azimuthal elastic impedance is given pri , Equation (8) can be solved by the least squares algorithm based on the model parameters and is expressed as:

[0099]

[0100] In the formula, represents the inverse of the model parameter covariance matrix, and \(\sigma\) represents the damping factor related to the signal-to-noise ratio. By solving using formula (10), the azimuthal elastic impedance data can be estimated from the azimuthal seismic data of partial angle stacking.

[0101] Step 2: After obtaining the azimuthal elastic impedance data, calculate the elastic impedance ratio between two different observation azimuths, and then use the Newton method to solve the normal weakness parameter and tangential weakness parameter of the fracture.

[0102] For two different observation azimuths (such as azimuth and azimuth ), the elastic impedance ratio can be expressed as:

[0103]

[0104] Among them,

[0105]

[0106] In formula (11), \(DEI\) represents the ratio of the elastic impedance data of two different azimuths.

[0107] Equation (11) can be expressed in matrix form as:

[0108] d f =G2m f (13)

[0109] Among them, \(G2\) represents the non-linear forward operator related to the incident angle and azimuth angle, and

[0110]

[0111] Similar to the mathematical and physical forward equation of Equation (13), its model parameters can be solved by the Newton method and expressed as:

[0112] m l+1 =m l +β l Δm l (15)

[0113] In the formula, β l represents the step size of each iteration and can be obtained by the line search method. m l+1 represents the predicted result of the model parameters after the l-th iteration, m l represents the initial value of the model parameters at the l-th iteration, Δm l respectively represent the perturbation amounts of the model parameters at the l-th iteration and can be expressed as:

[0114] Δm = -H -1 Y(16)

[0115] In the formula, Y and H respectively represent the first derivative and the second derivative of the model observation data vector d mod (m) with respect to the model parameter m and can be approximately solved as:

[0116]

[0117] In the formula, d input (m) represents the input actual observation data vector. The symbol "diag" represents a diagonal matrix.

[0118] In the traditional method, the values of Y and H are usually calculated by the finite difference method. In the present invention, they can be directly calculated using the derived elastic impedance equation. For the azimuthal elastic impedance difference equations (11) and (13) derived in the present invention, then in Equation (17) can be expressed as:

[0119]

[0120] where

[0121]

[0122] Substituting Equation (18) into Equations (15)-(17), the normal weakness parameter and the tangential weakness parameter of the fracture can be obtained.

[0123] In the third step, after obtaining the normal weakness parameter and the tangential weakness parameter of the fracture, using the derived ORT medium azimuthal elastic impedance equation, calculate the elastic impedance data related only to the incident angle, and also use the Newton method to invert properties A, property B, property C, and property D, and then calculate the anisotropy parameter ε using the inverted properties A and C.

[0124] As can be seen from Equation (5), after obtaining the normal weakness parameter and the tangential weakness parameter of the fracture, the elastic impedance equation related only to the incident angle can be expressed as:

[0125]

[0126] Equation (20) can be expressed in matrix form as:

[0127] d v = G3m v (21)

[0128] where G3 represents a non - linear forward operator that depends only on the incident angle, and

[0129]

[0130] According to equations (20) and (21), then in equation (17), can be expressed as:

[0131]

[0132] where,

[0133]

[0134] Substituting equation (23) into equations (15) - (17), reasonable predictions of the parameters of property A, property B, property C, and property D can be obtained. Then, the P - wave anisotropy parameter ε can be calculated using equation (7).

[0135] To test the feasibility and effectiveness of the inversion method proposed in this invention, synthetic tests were carried out using the time - domain logging curves of a certain fractured reservoir. Synthetic azimuthal seismic gather was generated by convolving the reflection coefficient equation (1) with a Ricker wavelet having a dominant frequency of 35 Hz, and Gaussian random noise with signal - to - noise ratios of 10:1 and 4:1 was added to the synthetic seismic gather respectively, as Figure 2 shown, where the incident angles of the synthetic seismic gather are 6°, 16°, 26°, and 36° respectively; the azimuth angles are 0° and 90°. Figure 3 and Figure 4 respectively show the azimuthal elastic impedance inversion results at signal - to - noise ratios of 10:1 and 4:1 ( Figure 3 and Figure 4 in which the solid line represents the true model, the dotted line represents the initial model, and the dashed line represents the inversion result), Figure 5 and Figure 6 respectively show the model parameter inversion results at signal - to - noise ratios of 10:1 and 4:1 ( Figure 5 and Figure 6 in which the solid line represents the true model, the dotted line represents the initial model, and the dashed line represents the inversion result). It can be observed from the inversion results that although the inversion results of the model parameters become worse as the signal - to - noise ratio decreases, overall, the inversion results of the model parameters still match well with the true values when the signal - to - noise ratio is 4:1. This shows that even in the presence of moderate noise, the step - by - step inversion method for orthotropic media based on the second - order derivative of elastic impedance proposed in this invention can still obtain a relatively reliable inversion result.

[0136] Field seismic data collected from a fractured reservoir work area in the western part of China are selected to further verify the feasibility and effectiveness of the inversion method proposed by the present invention. The lithology of the target layer in this work area is carbonaceous shale, and the lithologies of the upper and lower layers of the target layer are siltstone mudstone and limestone respectively. The target layer in the work area is a gas-bearing reservoir with low porosity and permeability. Geological core and logging data indicate that multiple vertical and nearly vertical microfractures are developed in the target layer. Therefore, the target reservoir can be equivalent to an ORT medium. The seismic azimuths used for inversion are 45°, 75°, and 135° respectively, and the incident angles corresponding to each azimuth are 8° (superposition of 4 - 12°), 16° (superposition of 12 - 20°), 24° (superposition of 20 - 28°), and 32° (superposition of 28 - 36°). The inversion method proposed by the present invention is used to estimate the model parameters, and the prediction results are as Figure 7 shown. The black curve in the figure represents the high-cut filtering display of the true model parameter curve at the well location. It can be seen that Property A (i.e., the longitudinal wave modulus), Property B (i.e., the anisotropic shear modulus), Property C (i.e., the anisotropic longitudinal wave modulus), and Property D (i.e., the density) are all at low values near the target layer section (about 1420 ms), while the normal fracture weakness, tangential fracture weakness, and anisotropy parameter are all at high values, which is basically consistent with the prior information of the work area. In addition, it can be observed from the figure that the inversion results have good continuity and they fit well with the well curves, which further confirms the feasibility and effectiveness of the inversion method proposed by the present invention.

[0137] The above-described embodiments are merely preferred embodiments given to fully illustrate the present invention, and the protection scope of the present invention is not limited thereto. Equivalent substitutions or transformations made by those skilled in the art based on the present invention are all within the protection scope of the present invention.

Claims

1. A step-by-step inversion method for orthotropic media based on the second-order derivative of elastic impedance, characterized in that: It includes the following steps: In the first step, the corresponding azimuthal elastic impedance data are estimated from the azimuthal seismic data stacked at certain angles; The second step is to calculate the ratio of the elastic impedance data at two different directions to solve the normal weakness parameter and tangential weakness parameter of the crack; The third step is to calculate the elastic impedance data containing only the incident angle information based on the parameter combination and mathematical approximation of the azimuthal elastic impedance equation of the orthotropic medium, invert the attributes A, B, C, and D, and use the ratio of attributes A to C to obtain the longitudinal wave anisotropy parameter ε. The azimuthal elastic impedance equation of the orthotropic medium based on parameter combination and mathematical approximation is: In formula (5), the symbol "exp" represents exponential operation, and the subscript symbol "0" represents the constant background medium elastic parameter. In equations (5) and (6), AEI represents the azimuthal elastic impedance, θ represents the incident angle, represents the azimuth angle, M and μ represent the longitudinal wave modulus and shear wave modulus of the isotropic background, ρ represents the density, g = μM; attribute A represents the longitudinal wave modulus; attribute B represents the anisotropic shear modulus. When the anisotropy parameter is 0, attribute B represents the isotropic background shear modulus; attribute C represents the anisotropic longitudinal wave modulus; attribute D represents the density; δ N Indicates the normal weakness parameter of the crack; δ T represents the tangential weakness parameter of the crack; ε represents the longitudinal wave anisotropy parameter.

2. The step-by-step inversion method for orthotropic media based on the second-order derivative of elastic impedance according to claim 1, characterized in that: In the first step: For seismic data containing i sampling points, the relationship between seismic data and azimuthal elastic impedance data is expressed in matrix form as follows: S PP =G1X(8), In formula (8), In formula (9), j and k represent the number of incident angles and azimuth angles, respectively, and s i represents the seismic record of the i-th sampling point, lnAEI represents the logarithm of the azimuthal elastic impedance, W represents the wavelet matrix, D o represents the difference operation matrix; When the initial model of azimuthal elastic impedance X is given pri Then, equation (8) is solved by the least squares algorithm based on model parameters and expressed as: In formula (10), represents the inverse of the covariance matrix of the model parameters, σ represents the damping factor related to the signal-to-noise ratio, and is solved using formula (10) to estimate the azimuthal elastic impedance data from the azimuthal seismic data stacked at certain angles.

3. The step-by-step inversion method for orthotropic media based on the second-order derivative of elastic impedance according to claim 1 or 2, characterized in that: In the second step: For two different observation positions and direction The elastic impedance ratio is expressed as: In formula (10), DEI represents the ratio of elastic impedance data in two different directions, Formula (11) is expressed in matrix form as: d f =G2m f (13), In formula (13), G2 represents the nonlinear forward operator related to the incident angle and azimuth angle, and The model parameters of formula (13) are solved by Newton method as follows: m l+1 =m l +b l Δm l (15), In formula (15), β l Indicates the step size of each iteration, obtained by line search method, m l+1 Represents the prediction results of the model parameters after the lth iteration, m l Indicates the initial value of the model parameters at the lth iteration, Δm l They represent the perturbation of the model parameters at the lth iteration as follows: Δm=-H -1 Y(16), In formula (16), Y and H represent the model observation data vector d mod (m) The first-order derivative and second-order derivative of the model parameter m are approximately solved as follows: In formula (16), d input (m) represents the actual observation data vector input, and the symbol "diag" represents a diagonal matrix; According to formulas (11) and (13), the formula (17) Expressed as: In formula (18), Substituting Equation (18) into Equations (15)-(17), the normal weakness parameter and tangential weakness parameter of the crack are obtained.

4. The step-by-step inversion method for orthotropic media based on the second-order derivative of elastic impedance according to claim 3, characterized in that: In the third step: The elastic impedance equation that depends only on the angle of incidence is expressed as: Formula (20) is expressed in matrix form as: d v =G3m v (21), In formula (20), G3 represents the nonlinear forward operator that is only related to the incident angle, and According to formulas (20) and (21), the formula (17) Expressed as: In formula (23), Substituting equation (23) into equations (15)-(17), we can obtain the parameters of attribute A, attribute B, attribute C, and attribute D; Utilization Calculate the longitudinal wave anisotropy parameter ε.

Images