A method for inversion of fracture weakness parameters

Abstract

The present invention discloses a fracture weakness parameter inversion method, comprising the following steps: S1: acquiring 3CVSP data of a target study area, and processing to obtain an upgoing S1 wave and an upgoing S2 wave; S2: obtaining the P wave velocity and the S wave velocity of the overlying stratum above the target layer according to the first arrival waves of the P wave and the S wave in the logging data and the 3CVSP data; S3: calculating and obtaining the travel time of the target layer at different detectors; S4: picking up the phase axis one closest to the target layer in the upgoing S1 wave according to the travel time, and obtaining the observed PS1 propagation time; S5: using the upgoing S1 wave, extracting the phase axis two of the target layer in the upgoing S2 wave through waveform matching, and obtaining the observed PS2 propagation time according to the travel time corresponding to the phase axis two; S6: establishing an objective function based on the observed PS1 propagation time and the observed PS2 propagation time, and inverting and obtaining the fracture weakness parameter according to the objective function. The present invention can accurately invert the fracture weakness parameter using the propagation time of the PS1 wave and the PS2 wave.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration and development, and in particular to a fracture weakness parameter inversion method. Background Art

[0002] The fractures and faults in fractured oil and gas reservoirs are of great significance to the storage and migration of oil and gas. Especially with the development of shale gas and ultra-deep oil and gas reservoirs, the development of fractures is becoming more and more important to the exploration and development of oil and gas reservoirs. When shear waves pass through fractured formations, they split into fast and slow shear waves. This effect is called shear wave splitting. The time difference and azimuth of shear wave splitting reflect the density and development of fractures in the underground medium. Through shear wave splitting analysis, relevant information about fractures can be obtained, thereby realizing the prediction of fracture development.

[0003] The response characteristics of seismic wave amplitudes to orthotropic (ORT) media can be used to invert fracture parameters. The propagation times of various waves in anisotropic media are usually an important tool for calculating fracture parameters. Many researchers have proposed different approximate methods for the propagation time of ORT media composed of a set of fractures. However, the quality of seismic data obtained from ground three-component observations, especially conversion waves, is difficult to guarantee the accuracy of amplitude and travel time inversion, which makes it difficult to apply the travel time approximation method based on ground seismic observation systems. Summary of the invention

[0004] In view of the above problems, the present invention aims to provide a method for inversion of fracture weakness parameters.

[0005] The technical solution of the present invention is as follows:

[0006] A fracture weakness parameter inversion method comprises the following steps:

[0007] S1: Acquire 3CVSP data of the target study area, and process the 3CVSP data to obtain upgoing S1 waves and upgoing S2 waves;

[0008] S2: Acquire well logging data of the target study area, and acquire P-wave velocity and S-wave velocity of the overlying strata above the target layer according to the well logging data and the first arrival waves of P-wave and S-wave in the 3CVSP data;

[0009] S3: Calculating the travel time of the target layer at different geophones according to the average velocity of the target layer in the logging data;

[0010] S4: picking up the event axis 1 closest to the target layer in the upgoing S1 wave according to the travel time, and obtaining the observed PS1 propagation time;

[0011] S5: using the upgoing S1 wave, extracting the second event of the target layer in the upgoing S2 wave through waveform matching, and obtaining the observed PS2 propagation time according to the travel time corresponding to the second event;

[0012] S6: Establishing an objective function based on the observed PS1 propagation time and the observed PS2 propagation time, and obtaining fracture weakness parameters by inverting the objective function.

[0013] Preferably, in step S6, the objective function is:

[0014]

[0015]

[0016] For the offset projection of PS1 wave on the X-axis and Y-axis:

[0017]

[0018] For the offset projection of PS2 wave on X and Y axis:

[0019]

[0020] Where: E is the objective function; X and Y are intermediate parameters; h is the depth of the detector; η1 and η2 are the weights of PS1 and PS2 waves respectively; T obs1 and T obs2 are the observed propagation time of PS1 and PS2 respectively; T' PS1 and T' PS2 The modeled PS1 propagation time and PS2 propagation time are respectively; the subscript m = PS1, PS2, representing the PS1 wave and PS2 wave respectively; C 0,m , C 20 ,m,C 02,m , C 40,m , C 22 ,m,C 04 ,m,N 20 ,m,N 02 ,m are coefficients; X P and Y P are the offset projections of the P wave on the X-axis and Y-axis respectively; S1 and Y S1 are the offset projections of the S1 wave on the X-axis and Y-axis respectively; S2 and Y S2 They are the offset projections of the S2 wave on the X-axis and Y-axis respectively.

[0021] Preferably, the offset projection of the P wave on the X-axis and the Y-axis is calculated by the following formula:

[0022]

[0023] Where: h0 is the thickness of the target layer; p 3,P is the vertical slowness of the P wave; p1 is the horizontal slowness in the x direction; Z is the depth of the bottom interface of the target layer; p 3,Pup is the vertical slowness of the P wave in the isotropic upper layer; p2 is the horizontal slowness in the y direction;

[0024] The offset projections of the S1 and S2 waves on the X and Y axes are calculated by the following formula:

[0025]

[0026] Where: X S and Y S are the offset projections of the S wave on the X-axis and Y-axis respectively; p 3,S is the vertical slowness of the S wave; p 3,Sup is the vertical slowness of S waves in the isotropic upper layer.

[0027] As a preference, the vertical slowness p 3,P and p 3,S Calculate using the following formula:

[0028] p 3,k =A 0,k +A 20,k p1 2 +A 02,k p2 2 +A 40,k p1 4 +A 04,k p2 4 +A 22,k p1 2 p2 2 (7)

[0029] Where: subscript k = S1, S2, P, representing S1 wave, S2 wave and P wave respectively; A 0,k , A 20,k , A 02,k , A 40,k , A 04,k , A 22,k All are coefficients;

[0030] The coefficient A is obtained by performing a 4th-order Taylor expansion as follows:

[0031]

[0032]

[0033] Where: P and Q are intermediate parameters; n = 0, 1, 2, representing S1 wave, S2 wave and P wave respectively; a, b, c, d are the elements c in the stiffness coefficient matrix ij Function expression of , where i, j = 1, 2…6;

[0034] The stiffness coefficient matrix is:

[0035]

[0036] The function expressions of a, b, c, and d are:

[0037]

[0038] Where: c 11b 、c 12b 、c 13b 、c 33b 、c 44b 、c 66b are the stiffness coefficients of background VTI; ΔN is the normal weakness; ρ is the density; ΔV is the vertical tangential weakness; ΔH is the horizontal tangential weakness.

[0039] Preferably, the vertical slowness of P and S waves in the isotropic upper layer is calculated by the following formula:

[0040]

[0041] Where: c 33up and c 44up are the squares of the P-wave and S-wave velocities, respectively.

[0042] Preferably, the coefficient C is obtained by the following steps:

[0043] First, we obtain the precise formulas for the travel time of the PS1 and PS2 waves, which are as follows:

[0044] T PS1 =h0p 3,P +(Z-h0)p 3,Pup +X P p1+Y P p2+h0p 3,S1 +(Z-h0-h)p 3,S1 +X S1 p1+Y S1 p2. (13)

[0045] T PS2 =h0p 3,P +(Z-h0)p 3,Pup +X P p1+Y P p2+h0p 3,S2+(Z-h0-h)p 3,S2 +X S2 p1+Y S2 p2. (14)

[0046] Where: T PS1 and T PS2 are the precise travel times of PS1 and PS2 waves respectively; h0 is the thickness of the target layer; p 3,P is the vertical slowness of the P wave; Z is the depth of the bottom interface of the target layer; p 3,Pup is the vertical slowness of the P wave in the isotropic upper layer; p1 is the horizontal slowness in the x direction; p2 is the horizontal slowness in the y direction; p 3,S1 is the vertical slowness of S1 wave; p 3,S2 is the vertical slowness of the S2 wave;

[0047] Then, T is obtained according to the precise formula of the travel time of PS1 wave and PS2 wave. m 2 The exact formula of

[0048] Finally, for the T m 2 The exact formula is expanded by the fourth order Taylor about p1 and p2, and the coefficients of the terms p1 zero p2 zero, p1 quadratic p2 zero, p1 zero p2 quadratic, p1 quartic p2 zero, p1 quadratic p2 quadratic, p1 zero p2 quartic are obtained as the C 0,m , C 20,m , C 02,m , C 40,m , C 22,m , C 04,m ;

[0049] The coefficient N is obtained by the following formula:

[0050]

[0051]

[0052] Where: M X,m and M Y,m All are intermediate parameters; X m , Y m are the offset projections of PS1 wave or PS2 wave on the X-axis and Y-axis respectively; T m It is the precise time of PS1 wave or PS2 wave.

[0053] The beneficial effects of the present invention are:

[0054] The present invention can accurately invert the fracture weakness parameters of ORT media from 3CVSP data; specifically, firstly, the data used is 3CVSP data, which has less interference with the upward shear wave, and the result of the shear wave travel time inversion observation will be more accurate compared with the ground multi-component observation; secondly, the present invention derives the approximate equations of the propagation time of PS1 and PS2 waves, and uses Taylor expansion to calculate the approximate coefficient of the vertical velocity; then, the rational form approximation is used to establish the travel time approximation about the shot point coordinates, which converts the calculation of the travel time from the complex phase velocity domain to the group velocity domain. Although the approximation assumes close to the offset, the numerical test results show that it has high accuracy for the offset close to the target layer depth; finally, the present invention establishes the objective function of inverting the fracture weakness parameters and anisotropy parameters based on the analytical travel time formula, and the actual data test shows that the fracture weakness parameters and anisotropy parameters have different directional sensitivities, which can be used to optimize the observation system in future research. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0056] Figure 1 is a schematic diagram of a VSP observation system according to a specific embodiment;

[0057] Figure 2 It is a schematic diagram of the ground shot point position of a VSP observation system of a specific embodiment; wherein the dots are the shot point positions;

[0058] Figure 3 It is a schematic diagram of the original VSP data of a specific embodiment; wherein the left figure is the R component and the right figure is the T component;

[0059] Figure 4 Schematic diagram of pure PS wave data after separation of original VSP data in a specific embodiment; wherein the left figure is PS1 wave data, and the right figure is PS2 wave data;

[0060] Figure 5 It is a schematic diagram of the inversion results of fracture weakness parameters of a specific embodiment. DETAILED DESCRIPTION

[0061] The present invention is further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, in the absence of conflict, the embodiments in this application and the technical features in the embodiments can be combined with each other. It should be noted that, unless otherwise specified, all technical and scientific terms used in this application have the same meanings as those generally understood by those of ordinary skill in the art to which this application belongs. The words "including" or "comprising" and the like used in the disclosure of the present invention mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects.

[0062] The present invention provides a method for inversion of fracture weakness parameters, comprising the following steps:

[0063] S1: Acquire 3CVSP data of the target study area, and process the 3CVSP data to obtain upgoing S1 waves and upgoing S2 waves.

[0064] S2: Acquire well logging data of the target study area, and acquire P-wave velocity and S-wave velocity of the overlying strata above the target layer based on the well logging data and the first arrival waves of P-wave and S-wave in the 3CVSP data.

[0065] S3: Calculate the travel time of the target layer at different detectors according to the average velocity of the target layer in the logging data.

[0066] S4: Pick up the event axis 1 closest to the target layer in the upgoing S1 wave according to the travel time to obtain the observed PS1 propagation time.

[0067] S5: Utilizing the upgoing S1 wave, extract the second event of the target layer in the upgoing S2 wave through waveform matching, and obtain the observed PS2 propagation time according to the travel time corresponding to the second event (the travel time corresponding to the second event is the observed PS2 propagation time).

[0068] S6: Establishing an objective function based on the observed PS1 propagation time and the observed PS2 propagation time, and obtaining fracture weakness parameters by inverting the objective function.

[0069] It should be noted that in steps S1-S5, separating the upgoing S1 wave and the upgoing S2 wave in the 3CVSP data, obtaining the wave velocity of the overlying stratum above the target layer according to the first arrival wave, calculating the travel time according to the average velocity, picking the event axis, waveform matching, etc. are all existing technologies, and the specific steps will not be repeated here.

[0070] In a specific embodiment, the objective function is:

[0071]

[0072]

[0073] For the offset projection of PS1 wave on the X-axis and Y-axis:

[0074]

[0075] For the offset projection of PS2 wave on X and Y axis:

[0076]

[0077] Where: E is the objective function; X and Y are intermediate parameters; h is the depth of the detector; η1 and η2 are the weights of PS1 and PS2 waves respectively; T obs1 and T obs2 are the observed propagation time of PS1 and PS2 respectively; T' PS1 and T' PS2 The modeled PS1 propagation time and PS2 propagation time are respectively; the subscript m = PS1, PS2, representing the PS1 wave and PS2 wave respectively; C 0,m , C 20,m , C 02,m , C 40,m , C 22,m , C 04,m 、N 20,m 、N 02,m All are coefficients; X P and Y P are the offset projections of the P wave on the X-axis and Y-axis respectively; S1 and Y S1 are the offset projections of the S1 wave on the X-axis and Y-axis respectively; S2 and Y S2 They are the offset projections of the S2 wave on the X-axis and Y-axis respectively.

[0078] It should be noted that, in the above embodiment, the X axis and the Y axis are the X axis and the Y axis in the rectangular coordinate system established by the geoid. The weights of the PS1 wave and the PS2 wave are determined according to the quality of the two signals. If the quality of the two signals is the same, both are 0.5.

[0079] In the above embodiment, the T' m The rational form approximate expression of the conversion wave travel time of the ground seismic observation system derived by the present invention is as follows:

[0080] A set of vertical cracks in the VTI background constitutes the ORT medium. After normalization by density, the stiffness coefficient matrix can be expressed as:

[0081]

[0082] Where: c 11b 、c12b 、c 13b 、c 33b 、c 44b 、c 66b are the stiffness coefficients of background VTI; ΔN is the normal weakness; ρ is the density; ΔV is the vertical tangential weakness; ΔH is the horizontal tangential weakness.

[0083] The anisotropy parameters of the background VTI medium are:

[0084]

[0085] Here the vertical velocities of P and S waves are:

[0086]

[0087] The slowness vector p = (p1, p2, p3) and displacement vector u of the ORT medium satisfy:

[0088] [Γ(p)-I]u=0 (19)

[0089] where I is the 3×3 unit diagonal matrix and Γ is the Christoffel matrix:

[0090]

[0091] If formula (19) has a solution, its determinant satisfies:

[0092] |Γ(p)-I|=0, (21)

[0093] Formula (21) can be written as vertical slowness p3 2 The cubic form of is:

[0094] ap3 6 +bp3 4 +cp3 2 +d=0 (22)

[0095] Where a, b, c, and d are the elements c in the stiffness coefficient matrix shown in formula (10) ij The function expressions are:

[0096]

[0097] p3 2 The analytical solution is:

[0098]

[0099]

[0100] Where: P and Q are intermediate parameters; n = 0, 1, 2, representing S1 wave, S2 wave and P wave respectively;

[0101] By using the fourth-order Taylor expansion of formula (8) and assuming a near zero offset, we obtain the vertical slowness approximation:

[0102] p 3,k =A 0,k +A 20,k p1 2 +A 02,k p2 2 +A 40,k p1 4 +A 04,k p2 4 +A 22,k p1 2 p2 2 (7)

[0103] Where: subscript k = S1, S2, P, representing S1 wave, S2 wave and P wave respectively; A 0,k , A 20,k , A 02,k , A 40,k , A 04,k , A 22,k All are coefficients;

[0104] In the VSP observation system, the converted shear wave originates from the downgoing P wave and the upgoing converted shear wave. The thickness of the target layer is h0, the depth of the bottom interface of the target layer is Z, and the depth of the detector is h. Based on this geometric shape, the present invention derives the corresponding kinematic parameters as follows:

[0105] The offset projection of the P wave on the X-axis and Y-axis is calculated by the following formula:

[0106]

[0107] Where: h0 is the thickness of the target layer; p 3,P is the vertical slowness of the P wave; p1 is the horizontal slowness in the x direction; Z is the depth of the bottom interface of the target layer; p 3,Pup is the vertical slowness of the P wave in the isotropic upper layer; p2 is the horizontal slowness in the y direction;

[0108] The offset projections of the S1 and S2 waves on the X and Y axes are calculated by the following formula:

[0109]

[0110] Where: X S and Y S are the offset projections of the S wave on the X-axis and Y-axis respectively; p 3,S is the vertical slowness of the S wave; p3,Sup is the vertical slowness of S waves in the isotropic upper layer.

[0111] In formulas (5)-(6), the vertical slowness of P waves and S waves in the isotropic upper layer satisfies:

[0112]

[0113] Where: c 33up and c 44up are the squares of the P-wave and S-wave velocities, respectively.

[0114] The exact travel time of the PS1 wave is:

[0115] T PS1 =h0p 3,P +(Z-h0)p 3,Pup +X P p1+Y P p2+h0p 3,S1 +(Z-h0-h)p 3,S1 +X S1 p1+Y S1 p2. (13)

[0116] The precise travel time of PS2 wave is:

[0117] T PS2 =h0p 3,P +(Z-h0)p 3,Pup +X P p1+Y P p2+h0p 3,S2 +(Z-h0-h)p 3,S2 +X S2 p1+Y S2 p2. (14)

[0118] Where: T PS1 and T PS2 are the precise travel times of PS1 and PS2 waves respectively; h0 is the thickness of the target layer; p 3,P is the vertical slowness of the P wave (obtained by Formula 7, i.e., k = P in Formula 7); Z is the depth of the bottom interface of the target layer; p 3,Pup is the vertical slowness of the P wave in the isotropic upper layer; p1 is the horizontal slowness in the x direction; p2 is the horizontal slowness in the y direction; p 3,S1 is the vertical slowness of the S1 wave (obtained by Formula 7, i.e., k = S1 in Formula 7); p 3,S2 is the vertical slowness of the S2 wave (obtained by Formula 7, i.e. k = S2 in Formula 7);

[0119] Formulas (13)-(14) represent the precise propagation time of the PS1 wave and the PS2 wave, and their expressions are functions of and , indicating that they are in the phase velocity domain. However, in practical applications, the coordinates of the shot point and the detector point are known and belong to the group velocity domain. Therefore, it is necessary to convert the expression in the phase velocity domain into the group velocity domain containing the shot point coordinates. Thus, the rational form approximate expression of the conversion wave travel time of the ground seismic observation system shown in formula (2) of the present invention is obtained.

[0120] In a specific embodiment, the coefficient C in formula (2) is obtained by the following steps:

[0121] First, obtain the precise formulas for the travel times of the PS1 and PS2 waves shown in equations (13)-(14);

[0122] Then, T is obtained according to the precise formula of the travel time of PS1 wave and PS2 wave. m 2 The exact formula of

[0123] Finally, for the T m 2 The exact formula is expanded by the fourth order Taylor about p1 and p2, and the coefficients of the terms p1 zero p2 zero, p1 quadratic p2 zero, p1 zero p2 quadratic, p1 quartic p2 zero, p1 quadratic p2 quadratic, p1 zero p2 quartic are obtained as the C 0,m , C 20,m , C 02,m , C 40,m , C 22,m , C 04,m ;

[0124] The coefficient N is obtained by the following formula:

[0125]

[0126]

[0127] Where: M X,m and M Y,m All are intermediate parameters; X m , Y m are the offset projections of PS1 wave or PS2 wave on the X-axis and Y-axis respectively; T m It is the precise time of PS1 wave or PS2 wave.

[0128] In a specific embodiment, the B area of ​​the T basin is taken as the target research area, and the fracture weakness parameter inversion method of the present invention is used to obtain its fracture weakness parameters. Figure 1 As shown in Figure 1. The distance between the shot point and the well is 1.8 km. The azimuth interval of each shot point relative to the well is 30°, as shown in Figure 1. Figure 2The depth range of the detectors is 3550-3750m, and the spacing is 10m. The original VSP shot detection set with an azimuth of 90° is shown in Figure 3 The raw data is a mixture of upgoing and downgoing P and S waves.

[0129] The target layer is a thin interbedded sandstone-mudstone sequence in the Triassic strata, with a depth range of 4050-4317 meters. The thin interbedded sandstone-mudstone structure constitutes the background VTI medium. The nearly vertical fractures developed therein form the ORT medium. The upper layer velocity is obtained by calculating the first arrival of the P wave and the S wave in the VSP record, which are 2.69km / S and 1.834km / S, respectively. The average velocities of the P wave and the S wave in the target layer are 3.64km / S and 1.91km / S, respectively. The upward PS1 and PS2 waves obtained by wave field separation are as follows: Figure 4 As shown in Figure 2, the main polarization direction of the fast shear wave obtained by shear wave splitting analysis is 35°. Therefore, Figure 2 The dashed arrows in the figure indicate the main azimuth of the crack, while the gray arrows indicate the direction perpendicular to the crack. These two axes (X and Y in the figure) constitute the two horizontal symmetry axes of the ORT medium.

[0130] To predict the weakness parameters of the cracks in the target layer, the velocity information is first used to calculate the travel time of the bottom of the target layer, e.g. Figure 4 The reflected wave of the target layer is shown by the left dashed line in the figure. Based on the similarity between fast shear waves and slow shear waves, the PS2 waves corresponding to these phase axes are identified and marked by dashed lines in the figure. Using the fracture weakness parameter inversion method of the present invention, according to Figure 2 The propagation time of PS waves at different shot points is used to calculate the crack and anisotropy parameters. The inversion results are shown in Figure 5 As shown, significant differences are shown for different azimuth angles.

[0131] Different parameters have different sensitivities at different azimuths. In order to obtain better results, the average results are calculated based on the sensitivity of these different parameters at different azimuths. First, according to the analysis of the symmetry axis direction (such as Figure 2 As shown), 125° and 305° are determined as the X-axis, and 35° and 215° are determined as the Y-axis. The sensitivity of ΔV in all azimuths is higher than that of other parameters, and the results at different azimuths are also relatively stable. Therefore, the average value of all directions is taken as the result. Since the errors of ΔN and ΔH are relatively small when the azimuth is close to 0° (X direction), the results of ΔN and ΔH near the X-axis (X±45°) are selected to calculate the average value. The results of the fracture weakness parameters are ΔV=0.153, ΔN=0.233, and ΔH=0.235. The average values ​​of the VTI parameter results near the Y-axis (±45°) are ε=0.144, γ=0.121, and δ=0.218.

[0132] In summary, the present invention can accurately invert and obtain the fracture weakness parameters of ORT media from 3CVSP data. Compared with the prior art, the present invention has significant progress.

[0133] The above description is only a preferred embodiment of the present invention and does not constitute any form of limitation to the present invention. Although the present invention has been disclosed as a preferred embodiment as above, it is not intended to limit the present invention. Any technician familiar with the profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the scope of the technical solution of the present invention.

Claims

1. A method for inversion of fracture weakness parameters, characterized in that: The following steps are involved: Step 1: Acquire 3CVSP data of the target study area, and process the 3CVSP data to obtain upgoing S1 waves and upgoing S2 waves; Step 2: Obtain the well logging data of the target study area, and obtain the P-wave velocity and S-wave velocity of the overlying strata above the target layer based on the well logging data and the first arrival waves of the P-wave and S-wave in the 3CVSP data; Step 3: Calculate the travel time of the target layer at different detectors according to the average velocity of the target layer in the logging data; Step 4: Pick the event axis 1 closest to the target layer in the upgoing S1 wave according to the travel time to obtain the observed PS1 propagation time; Step 5: using the upgoing S1 wave, extracting the second event of the target layer in the upgoing S2 wave through waveform matching, and obtaining the observed PS2 propagation time according to the travel time corresponding to the second event; Step 6: establishing an objective function based on the observed PS1 propagation time and the observed PS2 propagation time, and obtaining fracture weakness parameters by inverting the objective function; The objective function is: For the offset projection of PS1 wave on the X-axis and Y-axis: For the offset projection of PS2 wave on X and Y axis: Where: E is the objective function; X m and Y m is the intermediate parameter, and the subscripts m = PS1, PS2 represent the PS1 wave and PS2 wave respectively; when m = PS1, X PS1 and Y PS1 are the offset projections of the PS1 wave on the X-axis and Y-axis respectively; when m = PS2, X PS2 and Y PS2 are the offset projections of the PS2 wave on the X-axis and Y-axis respectively; h is the depth of the detector; η1 and η2 are the weights of the PS1 wave and the PS2 wave respectively; T obs1 and T obs2 are the observed propagation time of PS1 and PS2, respectively; T' PS1 and T' PS2 The modeled PS1 propagation time and PS2 propagation time are respectively; C 0,m , C 20,m , C 02,m , C 40,m , C 22,m , C 04,m 、N 20,m 、N 02,m All are coefficients; X P and Y P are the offset projections of the P wave on the X-axis and Y-axis respectively; S1 and Y S1 are the offset projections of the S1 wave on the X-axis and Y-axis respectively; S2 and Y S2 They are the offset projections of the S2 wave on the X-axis and Y-axis respectively.

2. The method for inversion of fracture weakness parameters according to claim 1, characterized in that: The offset projection of the P wave on the X-axis and Y-axis is calculated by the following formula: Where: h0 is the thickness of the target layer; p 3,P is the vertical slowness of the P wave; p1 is the horizontal slowness in the x direction; Z is the depth of the bottom interface of the target layer; p 3,Pup is the vertical slowness of the P wave in the isotropic upper layer; p2 is the horizontal slowness in the y direction; The offset projections of the S1 and S2 waves on the X and Y axes are calculated by the following formula: Where: X S and Y S are the offset projections of the S wave on the X-axis and Y-axis respectively; p 3,S is the vertical slowness of the S wave; p 3,Sup is the vertical slowness of S waves in the isotropic upper layer.

3. The method for inversion of fracture weakness parameters according to claim 2, characterized in that: Vertical slowness p 3,P and p 3,S Calculate using the following formula: p 3,k =A 0,k +A 20,k p1 2 +A 02,k p2 2 +A 40,k p1 4 +A 04,k p2 4 +A 22,k p1 2 p2 2 (7) Where: subscript k = S1, S2, P, respectively, represents S1 wave, S2 wave and P wave; when k = S1, p 3,S1 is the vertical slowness of the S1 wave; when k = S2, p 3,S2 is the vertical slowness of the S2 wave; A 0,k , A 20,k , A 02,k , A 40,k , A 04,k , A 22,k All are coefficients; The coefficient A is obtained by performing a 4th-order Taylor expansion as follows: Where: p3 is the vertical slowness; P and Q are intermediate parameters; n = 0, 1, 2, representing S1 wave, S2 wave and P wave respectively; a, b, c, d are the elements c in the stiffness coefficient matrix ij Function expression of, where i, j = 1, 2…6; i, j are the number of rows and columns of the stiffness coefficient matrix respectively; The stiffness coefficient matrix is: The function expressions of a, b, c, and d are: Where: c 11b 、c 12b 、c 13b 、c 33b 、c 44b 、c 66b are the stiffness coefficients of the background VTI; ΔN is the normal weakness; ρ is the density; ΔV is the vertical tangential weakness; ΔH is the horizontal tangential weakness.

4. The method for inversion of fracture weakness parameters according to claim 2, characterized in that: The vertical slowness of P and S waves in the isotropic upper layer is calculated by the following formula: Where: c 33up and c 44up are the squares of the P-wave and S-wave velocities, respectively.

5. The method for inversion of fracture weakness parameters according to any one of claims 1 to 4, characterized in that: The coefficient C is obtained by the following steps: First, we obtain the precise formulas for the travel time of the PS1 and PS2 waves, which are as follows: T PS1 =h0p 3,P +(Z-h0)p 3,Pup +X P p1+Y P p2+h0p 3,S1 +(Z-h0-h)p 3,S1 +X S1 p1+Y S1 p2. (13) T PS2 =h0p 3,P +(Z-h0)p 3,Pup +X P p1+Y P p2+h0p 3,S2 +(Z-h0-h)p 3,S2 +X S2 p1+Y S2 p2. (14) Where: T PS1 and T PS2 are the precise travel times of PS1 and PS2 waves respectively; h0 is the thickness of the target layer; p 3,P is the vertical slowness of the P wave; Z is the depth of the bottom interface of the target layer; p 3,Pup is the vertical slowness of the P wave in the isotropic upper layer; p1 is the horizontal slowness in the x direction; p2 is the horizontal slowness in the y direction; p 3,S1 is the vertical slowness of S1 wave; p 3,S2 is the vertical slowness of the S2 wave; Then, T is obtained according to the precise formula of the travel time of PS1 wave and PS2 wave. m 2 The exact formula of Finally, for the T m 2 The exact formula is expanded by the fourth order Taylor about p1 and p2, and the coefficients of the terms p1 zero p2 zero, p1 quadratic p2 zero, p1 zero p2 quadratic, p1 quartic p2 zero, p1 quadratic p2 quadratic, p1 zero p2 quartic are obtained as the C 0,m , C 20,m , C 02,m , C 40,m , C 22,m , C 04,m ; The coefficient N is obtained by the following formula: Where: M X,m and M Y,m All are intermediate parameters; X m , Y m are the offset projections of PS1 wave or PS2 wave on the X-axis and Y-axis respectively; T m It is the precise time of PS1 wave or PS2 wave.

Images