Inversion method based on three-dimensional reflection P wave monoclinic medium anisotropy parameters

Abstract

The invention discloses an inversion method based on three-dimensional reflection P wave monoclinic medium anisotropy parameters. The invention provides the three-dimensional reflection seismic wave travel time inversion method which is higher in stability, less in dependence on prior information and more sufficient in excavation of underground medium information. In the inversion process, firstly, dependence on prior information such as stratum thickness and vertical propagation speed is reduced, and meanwhile, more anisotropic parameters and relatively accurate stratum thickness can be obtained through inversion. Therefore, three-dimensional reflection seismic wave data can be explored more fully, richer underground medium information is provided, and a foundation is laid for a subsequent oil-gas exploration process, especially an unconventional oil-gas exploration process. A geophysical data processor can process the multi-azimuth single-layer horizontal plane reflection P wave travel time according to the method, anisotropy information of a bottom layer is obtained, or a multi-layer reflection P wave travel time and converted wave travel time inversion method is further developed according to the method provided by the invention.

Description

Technical Field

[0001] The present application relates to the field of petroleum geophysical exploration, and in particular to an inversion method for anisotropic parameters of monoclinic media based on three-dimensional reflection P waves. Background Art

[0002] Anisotropy is widespread in subsurface media due to factors such as the oriented arrangement of crystalline minerals, fractures, and the influence of geostress. In seismic exploration, this anisotropy manifests itself primarily as directional variations in seismic wave travel time, amplitude, and attenuation. The travel time of horizontal interface reflection waves is a crucial component of seismic response, and its forward modeling and inversion methods have been extensively studied.

[0003] In a single-azimuth observation system that ignores horizontal azimuth, reflected wave traveltimes vary only with offset, and their anisotropy manifests as a non-hyperbolic variation in the traveltime curve. The anisotropy parameters involved can only reflect the anisotropy of seismic wave propagation within a single vertical profile. Therefore, anisotropy parameters derived from single-azimuth traveltime inversion can effectively restore the anisotropy only when the medium is a transversely isotropic (VTI) medium with a vertical axis of symmetry. In a three-dimensional observation system containing multi-azimuth data, the addition of azimuthal anisotropy to reflected wave traveltimes—that is, the reflected wave traveltimes vary not only with offset but also with azimuth—supports the inversion of anisotropic media with azimuthal anisotropy, such as horizontally isotropic (HTI) media with a horizontal axis of symmetry and orthorhombic anisotropy. The development of high-density seismic technology and pre-processing technologies represented by OVT (Offset Vector Tile) technology can describe the azimuth-varying properties of monoclinic media, which have less symmetry and require more anisotropic parameters to describe, through reflection wave travel time inversion.

[0004] Research on monoclinic media began in 2000. Building on research on orthotropic anisotropic media, Grechka defined anisotropic media equivalent to multiple sets of vertical fractures as monoclinic media. He then inverted the anisotropic parameters using the NMO velocities of P- and S-waves reflected from their horizontal interfaces. Subsequent approximations and inversion methods for monoclinic traveltimes were largely based on this approach. While the inversion of NMO velocities can accurately determine the azimuthal properties of anisotropic media, it lacks the ability to describe properties such as anisotropy intensity. Furthermore, the NMO velocity inversion process relies heavily on prior information, such as the need for an axis or plane of the global coordinate system to coincide with a symmetry axis or plane of the medium, the thickness of the reflection layer, and the vertical P- and S-wave propagation velocities. Consequently, this inversion process has certain limitations. Summary of the Invention

[0005] In order to solve the above technical problems or at least partially solve the above technical problems, the present application provides an inversion method and device based on the anisotropic parameters of a three-dimensional reflected P-wave monoclinic medium.

[0006] According to one aspect of an embodiment of the present application, a method for inversion based on three-dimensional reflected P-wave anisotropy parameters is provided, comprising:

[0007] Determining a spatial coordinate system for the seismic data, wherein the spatial coordinate system includes establishing a three-dimensional global coordinate system and a two-dimensional local coordinate system;

[0008] Inverting the anisotropy parameter in the two-dimensional local coordinate system using the two-dimensional travel time data to obtain the first anisotropy parameter of the monoclinic medium;

[0009] Based on a pre-established conversion relationship of anisotropy parameters between a three-dimensional global coordinate system and a two-dimensional local coordinate system, the first anisotropy parameter is converted into a second anisotropy parameter in the three-dimensional global coordinate system, and the global anisotropy parameter is inverted based on the second anisotropy parameter to obtain the second anisotropy parameter;

[0010] A scanning calculation of the formation thickness is performed according to the second anisotropy parameter to obtain a calculation result.

[0011] According to another aspect of an embodiment of the present application, an inversion device based on three-dimensional reflected P-wave anisotropy parameters is provided, comprising:

[0012] A determination module, configured to determine a spatial coordinate system of seismic data, wherein the spatial coordinate system includes establishing a three-dimensional global coordinate system and a two-dimensional local coordinate system;

[0013] A local inversion module is used to invert anisotropic parameters in the two-dimensional local coordinate system using the two-dimensional travel time data to obtain a first anisotropic parameter of the monoclinic medium;

[0014] a global inversion module, configured to convert the first anisotropy parameter into a second anisotropy parameter in the three-dimensional global coordinate system based on a pre-established conversion relationship of anisotropy parameters between the three-dimensional global coordinate system and the two-dimensional local coordinate system, and perform an inversion of the global anisotropy parameter based on the second anisotropy parameter to obtain the second anisotropy parameter;

[0015] The calculation module is used to scan and calculate the formation thickness according to the second anisotropy parameter to obtain the target formation thickness of the underground medium.

[0016] According to another aspect of an embodiment of the present application, a storage medium is further provided, which includes a stored program, and the above steps are executed when the program is run.

[0017] According to another aspect of an embodiment of the present application, an electronic device is also provided, including a processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory communicate with each other through the communication bus; wherein: the memory is used to store computer programs; the processor is used to execute the steps in the above method by running the program stored in the memory.

[0018] An embodiment of the present application also provides a computer program product comprising instructions, which, when executed on a computer, enables the computer to execute the steps in the above method.

[0019] The above technical solution provided by the embodiment of the present application has the following advantages compared with the prior art:

[0020] The present application provides a three-dimensional reflection seismic wave travel time inversion method that is more stable, less dependent on prior information, and more comprehensive in mining underground medium information. In the inversion process, first of all, the dependence on prior information such as stratum thickness and vertical propagation velocity is reduced, and at the same time, more anisotropy parameters and relatively accurate stratum thickness can be obtained through inversion. This can more fully explore three-dimensional reflection seismic wave data, provide richer underground medium information, and lay the foundation for subsequent oil and gas exploration processes, especially unconventional oil and gas exploration processes. Geophysical data processing personnel can process multi-directional single-layer horizontal plane reflection P-wave travel time according to this method to obtain the anisotropy information of the bottom layer, or further develop multi-layer reflection P-wave travel time and conversion wave travel time inversion methods based on the method provided by the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present application and, together with the description, serve to explain the principles of the present application.

[0022] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, for ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0023] Figure 1 A flowchart of an inversion method based on three-dimensional reflected P-wave anisotropy parameters provided in an embodiment of the present application;

[0024] Figure 2 Schematic diagram of the reflected ray and global and local coordinate systems defined in the present invention;

[0025] Figure 3 A distribution range diagram of the seismic survey lines in Example 1 of the present application;

[0026] Figure 4 The covariance matrix for evaluating the accuracy of the anisotropic parameters in Example 1 of the present application;

[0027] Figure 5 A comparison chart of the anisotropy parameters obtained by inversion in Example 1 of the present application and the exact values;

[0028] Figure 6 This is a graph showing the relative travel time error as a function of the thickness of the preset reflective layer in Example 1 of the present application.

[0029] Figure 7 The covariance matrix for evaluating the accuracy of the anisotropic parameters in Example 2 of the present application;

[0030] Figure 8 A comparison diagram of the anisotropy parameters obtained by inversion in Example 2 of the present application and the exact values;

[0031] Figure 9 This is a graph showing the change in relative travel time error with the thickness of the preset reflective layer in Example 2 of this application.

[0032] Figure 10 A distribution range diagram of the seismic survey lines in Example 3 of the present application;

[0033] Figure 11 The covariance matrix for evaluating the accuracy of the anisotropic parameters in Example 3 of the present application;

[0034] Figure 12 A comparison diagram of the anisotropy parameters obtained by inversion in Example 3 of the present application and the exact values;

[0035] Figure 13 is a graph showing the change in relative travel time error with the thickness of the preset reflective layer in Example 3 of the present application;

[0036] Figure 14 A distribution range diagram of the seismic survey lines in Example 4 of the present application;

[0037] Figure 15 The covariance matrix for evaluating the accuracy of the anisotropic parameters in Example 4 of the present application;

[0038] Figure 16 A comparison diagram of the anisotropy parameters obtained by inversion in Example 4 of the present application and the exact values;

[0039] Figure 17 This is a graph showing the change in relative travel time error with the thickness of the preset reflective layer in Example 4 of this application.

[0040] Figure 18 A distribution range diagram of the seismic survey lines in Example 5 of the present application;

[0041] Figure 19 is the covariance matrix for evaluating the accuracy of the anisotropic parameters in Example 5 of the present application;

[0042] Figure 20 A comparison diagram of the anisotropy parameters obtained by inversion in Example 5 of the present application and the exact values;

[0043] Figure 21 is a graph showing the change in relative travel time error with the thickness of the preset reflective layer in Example 5 of the present application;

[0044] Figure 22 A schematic diagram of an inversion device based on three-dimensional reflected P-wave anisotropy parameters provided in an embodiment of the present application. DETAILED DESCRIPTION

[0045] In order to make the purpose, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all of the embodiments. The schematic embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation of the present application. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application.

[0046] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another similar entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, article, or device comprising a series of elements includes not only those elements, but also other elements not explicitly listed, or elements inherent to such process, method, article, or device. In the absence of further limitations, an element defined by the phrase "comprising a ..." does not exclude the presence of other identical elements in the process, method, article, or device comprising the element.

[0047] The present invention provides an inversion method for anisotropic parameters of a monoclinic medium based on three-dimensional P-wave reflection. The method provided by the present invention can be applied to any desired electronic device, such as a server or terminal. For ease of description, the method is not specifically limited here and is referred to as an electronic device.

[0048] According to one aspect of an embodiment of the present application, a method embodiment of an inversion method based on three-dimensional reflected P-wave anisotropy parameters is provided. Figure 1 A flowchart of an inversion method based on three-dimensional reflected P-wave anisotropy parameters provided in an embodiment of the present application is shown in FIG. Figure 1 As shown, the method includes:

[0049] Step S11 , determining the spatial coordinate system of the seismic data, wherein the spatial coordinate system includes establishing a three-dimensional global coordinate system and a two-dimensional local coordinate system.

[0050] In the embodiments of this application, Figure 2 As shown, first set up a cylindrical coordinate system is a three-dimensional global coordinate system, where the z-axis is vertically downward, The polar coordinate system formed by the above equations constitutes an axial plane perpendicular to the z-axis. The background medium is a stratum composed of a uniform monoclinic medium with a thickness of H. The z-axis of the three-dimensional global coordinate system is parallel to the intersection line of the fracture. The medium is centrally symmetrical with respect to the z-axis. The axial plane is Perpendicular to the crack surface, the medium is referenced to the axial plane The horizontal reflection plane at the bottom of the formation is symmetrical to the axial plane. In this case, the incident downgoing and reflected upgoing rays of the P wave lie in planes perpendicular to the reflecting surface and are symmetrical with respect to the normal to the reflecting surface.

[0051] According to the seismic wave observation system, a two-dimensional local coordinate system is defined. The seismic wave is excited from the source located at the origin S of the three-dimensional global coordinate system, reflected by the horizontal reflection surface, and received by the detector located at R. The origin of the two-dimensional local coordinate system is consistent with the source S. The straight line where the source point S and the detector R are located is defined as the x'1 axis of the two-dimensional local coordinate system, and the z axis of the three-dimensional global coordinate system is defined as the x'3 axis of the two-dimensional local coordinate system. The azimuth angle between the plane of the two-dimensional local coordinate system (x1', x3') and the three-dimensional global coordinate system is The reflected wave ray is parallel to the vertical plane of the 2D local coordinate system (x1', x3').

[0052] Step S12: Invert the anisotropic parameters in a two-dimensional local coordinate system using the two-dimensional travel time data to obtain the first anisotropic parameters of the monoclinic medium.

[0053] In the embodiment of the present application, an inversion of anisotropy parameters is performed in a two-dimensional local coordinate system using two-dimensional travel time data to obtain the first anisotropy parameter of the monoclinic medium, including the following steps A1-A3:

[0054] Step A1: construct an approximate formula for the variation of reflected P-wave travel time with offset distance in a two-dimensional local coordinate system.

[0055] In the embodiment of the present application, in the local coordinate system (x1', x3') plane, the reflected wave travel time T can be calculated by the following formula:

[0056]

[0057] Where x is the offset (the distance between the source S and the receiver R), H is the depth of the horizontal reflection plane (coinciding with the symmetry plane of the formation), T = T(x) represents the travel time of the reflected P wave, v = v(n) represents the ray velocity, and n is the phase vector.

[0058] Then, the reflection travel time formula is preprocessed. First, the expression method widely used in reflection wave time difference processing is used to simplify the above formula:

[0059]

[0060] in, is the offset distance normalized by the depth of the reflecting surface, α is the P-wave phase velocity in the reference isotropic medium, and T0 is the zero-offset two-way travel time in the reference isotropic medium. Formula (1) can be rewritten as:

[0061]

[0062] In formula (3), the ray velocity v(n) of anisotropic media is difficult to obtain, so the weak anisotropy approximation is introduced, and the ray velocity is approximately expressed by the phase velocity of the medium. The square of the phase velocity is used, that is, v 2 (N)~c 2 (N) approximates the ray velocity v(N). In the local coordinate system, the azimuth angle is not considered. The group velocity of a monoclinic medium can be expressed using the anisotropy parameter as:

[0063]

[0064] In formula (4), N = (N1, N2, N3) is the direction vector of the phase velocity. In the local coordinate system, N can be obtained by express""

[0065]

[0066] The anisotropic parameters in formula (4) include εx , ε z , η y , where ε x , ε z Respectively represent the difference between the P-wave phase velocity of the medium in the coordinate axis x1 and x3 directions and the reference velocity α; η y Represents the gradient of the P-wave velocity in the medium in the coordinate plane (x1, x3). The anisotropy parameter in the local coordinate system is related to the stiffness coefficient in the local coordinate system:

[0067]

[0068] Among them, A ij is the element in the density-normalized stiffness coefficient matrix of the formation medium in the local coordinate system, α is the P-wave phase velocity in the reference isotropic medium, and the value of α can be set arbitrarily without affecting the value of the square of the phase velocity in formula (4).

[0069] Substituting formula (4) into formula (3), we can get the following formula:

[0070]

[0071]

[0072] Where x is the offset, T = T0(x) represents the travel time of the reflected P wave, and the anisotropic parameters include ε x '、ε z '、η y ',ε x '、ε z 'represents the difference between the P-wave phase velocity of the medium in the coordinate axis x1 and x3 directions and the reference velocity α; η y 'Represents the gradient of the P-wave velocity change in the medium in the coordinate plane (x1, x3).

[0073] Step A2: construct the first inversion determinant using the approximate formula for the variation of reflected P-wave travel time with offset distance.

[0074] In the embodiment of the present application, the inversion determinant is constructed

[0075] By moving the known information to one side of the equal sign, formula (7) and formula (8) can be rewritten as:

[0076]

[0077] According to formula (9), the first inversion determinant is constructed as follows:

[0078] Gm=d, (10)

[0079] Where m=(ε x ',ηy ',ε x ') T , G is an N×M coefficient matrix, N is the number of detectors, M is the number of A parameters involved, d is a vector, and the row vector in the matrix G is of the following form:

[0080]

[0081] m consists of three parameters and has the following form:

[0082] m=(ε x , η y , ε z ) T , (12)

[0083] The vector d has the following form:

[0084]

[0085] is the normalized offset of the i-th detection point.

[0086] Step A3: Solve the first inversion determinant using the least squares method to obtain the first anisotropy parameter of the monoclinic medium.

[0087] Formula (10) can be solved by the least squares method, and the inversion objective function is:

[0088] ||Gm-d||2<κ, (13)

[0089] Where κ is the parameter that controls the inversion accuracy. Formula (13) is solved by formula (14):

[0090]

[0091] in, is the inversion operator, Among them, G T is the transpose of the matrix G, I is the N×N identity matrix, and ε is the stability factor. The larger ε is, the more stable the inversion process is, but the accuracy is lower.

[0092] Step S13: Based on the pre-established conversion relationship of anisotropic parameters between the three-dimensional global coordinate system and the two-dimensional local coordinate system, the first anisotropic parameter is converted into a second anisotropic parameter in the three-dimensional global coordinate system, and the global anisotropic parameter is inverted based on the second anisotropic parameter to obtain the second anisotropic parameter.

[0093] In the embodiments of this application,

[0094] In the global coordinate system, the azimuth angle needs to be considered The three-dimensional travel time comes from the seismic lines distributed in different directions on the plane, and the strike angle of the line is A strike angle is The anisotropic parameters of the local coordinate system composed of the survey lines can be converted into the anisotropic parameters of the global coordinate system through Bond transformation. The transformation matrix R in Bond transformation is:

[0095]

[0096] The anisotropic parameters of monoclinic media in the local coordinate system are represented by the superscript “′”, and the anisotropic parameters of orthotropic media in the global coordinate system are distinguished.

[0097] According to the 6×6 stiffness coefficient matrix conversion rule in Bond transformation, the anisotropy parameter ε of the orthotropic medium in the local coordinate system is x ',η y ',ε x 'It can be converted to the global coordinate system by the following formula:

[0098]

[0099]

[0100] ε′ z =ε z ; (19)

[0101] The second anisotropy parameter (global anisotropy parameter) includes ε x ,ε y ,ε z ,η x ,η y ,η z , where ε x ,ε y ,ε z Respectively represent the difference between the P-wave phase velocity of the medium in the coordinate axes x1, x2, and x3 and the reference velocity α; η x ,η y ,η z represent the change gradient of the P-wave velocity in the coordinate planes (x2, x3), (x1, x3), and (x1, x2), respectively; χ2 represents the degree to which the P-wave deviates from the axial direction in the plane (x1, x2); ξ 16 ,ξ 26 Represents the degree of deviation between the P wave direction and the fast and slow S wave directions in the plane (x1, x2).

[0102] In the embodiment of the present application, the anisotropy parameter in the global coordinate system is related to the stiffness coefficient in the global coordinate system:

[0103]

[0104] Among them, A ij is the element in the density-normalized stiffness coefficient matrix of the formation medium in the global coordinate system, and α is the P-wave phase velocity in the reference isotropic medium.

[0105] In the embodiment of the present application, inversion of the global anisotropy parameter is performed based on the second anisotropy parameter to obtain the second anisotropy parameter, including the following steps B1-B2:

[0106] Step B1: constructing a second inversion determinant based on the conversion relationship of anisotropy parameters between the three-dimensional global coordinate system and the two-dimensional local coordinate system.

[0107] In the embodiment of the present application, the second inversion determinant is as follows:

[0108] G1m1=d1;

[0109] G2m2=d2;

[0110] Among them, G1 is an N×6 coefficient matrix, N is the number of seismic lines distributed along different directions, G2 is an N×5 coefficient matrix, the number of rows N is the number of azimuths of the local coordinate system involved in the inversion, and the forms of vectors m1 and m2 are as follows:

[0111] m1=(ε x , ε y , η z , χ z ,ξ 16 ,ξ 26 ) T ;

[0112] m2=(η x , η y , η z ,ξ 16 ,ξ 26 ) T ;

[0113] The forms of vectors d1 and d2 are as follows:

[0114]

[0115]

[0116] In order to make the inverse problem of Equations (21) and (22) well-posed, the coefficient matrices G1 and G2 must be positive definite or overdetermined. Therefore, to obtain the A parameters in the global coordinate system by inverting the A parameters in the local coordinate system of multiple azimuth seismic lines, seismic reflection wave travel time data in at least six azimuths are required.

[0117] Step B2: Solve the second inversion determinant using the least squares method to obtain the second anisotropy parameter.

[0118] Formulas (24) and (25) can be solved by the least squares method, and the inversion objective function is:

[0119] ||G1m1-d1||2<κ, (29)

[0120] ||G2m2-d2||2<κ, (30)

[0121] Where κ is a parameter that controls the inversion accuracy. Formula (29) and (30) are solved by formula (31):

[0122]

[0123] in, is the inversion operator, i=1, 2 is the serial number of the inversion determinant.

[0124]

[0125] Among them, G i T is the matrix G i is the transpose of , I is the N×N identity matrix, ε is the stability factor, and the larger ε is, the more stable the inversion process is, but the accuracy is lower.

[0126] Step S14: performing a scanning calculation of the formation thickness according to the second anisotropy parameter to obtain a calculation result.

[0127] In an embodiment of the present application, after inverting the global anisotropy parameter based on the second anisotropy parameter to obtain the second anisotropy parameter, the method further includes: calculating the residual between the model synthetic data and the actual seismic data during the inversion process; determining the inversion accuracy of the inversion of the global anisotropy parameter by the residual; and when the reaction accuracy meets the preset conditions, performing a scanning calculation of the stratum thickness based on the second anisotropy parameter and the seismic data to obtain spatial information of the underground medium.

[0128] Specifically, to evaluate the accuracy of the inversion, it is necessary to calculate the model covariance matrix Gm of the inversion parameters through the following steps. First, the residual parameters are calculated using the residual statistics method:

[0129]

[0130] in: The matrix G is the inversion determinant G in step 3. The specific format is shown in formula (11). is obtained by calculating the inversion results in step 3 and step 5, and v is the difference between the number of rows and columns of matrix G.

[0131] The data covariance matrix is ​​then calculated based on the residual parameters:

[0132] C d ~σ 2 I, (34)

[0133] Where I is a 9×9 identity matrix. Finally, the data covariance matrix is ​​converted into the model covariance matrix:

[0134] C m =HC d H T , (35)

[0135] Here The model covariance matrix Gm represents the random distribution probability of the desired model parameter vector m. If randomness does not exist, the covariance matrix is ​​0, and the model parameter vector m has no error. The model covariance matrix Gm is a 9×9 square matrix, where the diagonal elements are the standard deviations of the corresponding anisotropic parameters, and the off-diagonal elements represent the correlation between the two related anisotropic parameters.

[0136] In an embodiment of the present application, a scanning calculation of the stratum thickness is performed based on the second anisotropy parameter and seismic data to obtain spatial information of the underground medium, including the following steps: obtaining a preset stratum thickness; calculating the first travel time of each detection point based on the second anisotropy parameter and the preset stratum thickness; obtaining the true second travel time of each detection point during the inversion process, and calculating the average relative travel time based on the first travel time and the second travel time; and selecting the stratum thickness with the smallest average relative travel time error as the calculation result.

[0137] Calculate the average relative travel time error T during a single scan AV (H), calculated as follows:

[0138]

[0139] The method provided in the embodiments of the present application establishes a three-dimensional global coordinate system and a two-dimensional local coordinate system, allowing for a more accurate description of the spatial location of seismic data. This helps improve the accuracy of data processing and analysis. By inverting the anisotropy parameters using two-dimensional travel-time data in a two-dimensional local coordinate system, the first anisotropy parameter of the monoclinic medium can be obtained. This helps obtain more detailed and accurate information about the subsurface medium. By using the pre-established conversion relationship between the anisotropy parameters of the three-dimensional global coordinate system and the two-dimensional local coordinate system, the first anisotropy parameter can be converted into a second anisotropy parameter in the three-dimensional global coordinate system. This allows for better integration of data from different coordinate systems and comprehensive analysis. Based on the second anisotropy parameter, the global anisotropy parameter can be inverted. This helps obtain more comprehensive and detailed anisotropy parameter information about the subsurface medium, improving the accuracy of seismic analysis. Based on the second anisotropy parameter, stratum thickness can be scanned and calculated. This helps obtain information about the thickness of subsurface strata, providing basic data and reference for applications such as geological exploration and earthquake prediction.

[0140] Based on this, the present application provides a three-dimensional reflection seismic wave travel time inversion method that is more stable, less dependent on prior information, and more comprehensive in mining underground medium information. In the inversion process, first of all, the dependence on prior information such as stratum thickness and vertical propagation velocity is reduced, and at the same time, more anisotropy parameters and relatively accurate stratum thickness can be obtained through inversion. This can more fully explore three-dimensional reflection seismic wave data, provide richer underground medium information, and lay the foundation for subsequent oil and gas exploration processes, especially unconventional oil and gas exploration processes. Geophysical data processing personnel can process multi-directional single-layer horizontal plane reflection P-wave travel time according to this method to obtain the anisotropy information of the bottom layer, or further develop multi-layer reflection P-wave travel time and conversion wave travel time inversion methods based on the method provided by the present invention.

[0141] Example 1

[0142] In a specific embodiment of the present invention, Table 1 lists the elastic parameters of the monoclinic medium used in the implementation of the present invention. These elastic parameters comprise a 6×6 stiffness coefficient matrix, which is symmetric along the diagonal elements. Therefore, only the upper triangular matrix is ​​shown. In this embodiment, the reflector layer thickness is 0.7 km. The travel times of the horizontal reflection waves in the three-dimensional medium used in this embodiment were calculated using the ANRAY program package and used as the actual measured seismic data.

[0143]

[0144] Table 1 Model elastic parameters used in Example 1

[0145] Table 2 shows the reference P-wave velocity used in Example 1 of the present invention and the precise value of the anisotropy parameter calculated according to formula (20) based on the reference P-wave velocity and the elastic parameters in Table 1. The reference velocity can be defined as an arbitrary value. Changes in the reference velocity will only affect the anisotropy parameter itself and will not affect the accuracy of the longitudinal wave velocity described by the anisotropy parameter and the corresponding reflection process. Table 2 also shows the distribution range of the observation system in the present invention, where is the strike angle of the seismic line, and x is the length of the seismic line.

[0146]

[0147] Table 2 Observation system range and exact values ​​of anisotropy parameters in Example 1

[0148] Figure 3 The distribution range of the seismic survey line in Example 1 is drawn according to the data distribution range parameters in Table 2. Figure 4 is the covariance matrix calculated according to formula (33)-formula (35). The diagonal elements of the matrix represent the variance of the corresponding parameters, and the non-diagonal elements represent the correlation between the corresponding two parameters. Figure 5 The comparison of the anisotropy parameters calculated by inversion in Example 1 with the exact values ​​in Table 2, where the error bars are Figure 4 The square roots of the diagonal elements of the covariance matrix in . Figure 6 This is a graph showing the average relative travel time error in step 7 as a function of the preset H.

[0149] Example 2

[0150] In a specific embodiment of the present invention, Table 3 shows the elastic parameters of the monoclinic medium used in the implementation of the present invention. These elastic parameters comprise a 6×6 stiffness coefficient matrix, which is symmetric along the diagonal elements. Therefore, only the upper triangular matrix is ​​shown. In this embodiment, the reflector thickness is 0.7 km. The travel times of the horizontal reflection waves in the three-dimensional medium used in this embodiment were calculated using the ANRAY program package. After adding 6 ms of random Gaussian noise, the data was used as the actual measured seismic data.

[0151]

[0152] Table 3 Model elastic parameters used in Example 2

[0153] Table 3 shows the reference P-wave velocity used in Example 2 of the present invention and the precise value of the anisotropy parameter calculated according to formula (20) based on the reference P-wave velocity and the elastic parameters in Table 3. The reference velocity can be defined as an arbitrary value. Changes in the reference velocity will only affect the anisotropy parameter itself and will not affect the accuracy of the longitudinal wave velocity described by the anisotropy parameter and the corresponding reflection process. Table 4 also shows the distribution range of the observation system in the present invention, where is the strike angle of the seismic line, and x is the length of the seismic line.

[0154]

[0155] Table 4 Observation system range and anisotropy parameter exact values ​​in Example 2

[0156] The distribution range of the seismic survey lines in this embodiment is the same as that in embodiment 1. Figure 3 , Figure 7 is the covariance matrix inverted in this embodiment calculated according to formula (33)-formula (35), the diagonal elements of the matrix represent the variance of the corresponding parameters, and the non-diagonal elements represent the correlation between the corresponding two parameters. Figure 8 The anisotropy parameters calculated by inversion in Example 2 are compared with the exact values ​​in Table 4, where the error bars are Figure 7 The square roots of the diagonal elements of the covariance matrix in . Figure 9 This is a graph showing the average relative travel time error in step 7 as a function of the preset H.

[0157] Example 3

[0158] In a specific embodiment of the present invention, Table 5 shows the elastic parameters of the monoclinic medium used in the implementation of the present invention. These elastic parameters comprise a 6×6 stiffness coefficient matrix, which is symmetric along the diagonal elements. Therefore, only the upper triangular matrix is ​​shown. In this embodiment, the reflector layer thickness is 0.7 km. The travel times of the horizontal reflection waves in the three-dimensional medium used in this embodiment were calculated using the ANRAY program package and used as the actual measured seismic data.

[0159] Table 5 Elastic parameters of the model used in Example 3

[0160]

[0161] Table 6 shows the reference P-wave velocity used in Example 3 of the present invention and the precise value of the anisotropy parameter calculated according to formula (20) based on the reference P-wave velocity and the elastic parameters in Table 5. The reference velocity can be defined as an arbitrary value. Changes in the reference velocity will only affect the anisotropy parameter itself and will not affect the accuracy of the longitudinal wave velocity described by the anisotropy parameter and the corresponding reflection process. Table 6 also shows the distribution range of the observation system in the present invention, where is the strike angle of the seismic line, and x is the length of the seismic line.

[0162]

[0163] Table 6 Observation system range and exact values ​​of anisotropy parameters in Example 3

[0164] Figure 10 The distribution range of the seismic survey lines in Example 3 is drawn according to the data distribution range parameters in Table 6. Figure 11 is the covariance matrix calculated according to formula (33)-formula (35). The diagonal elements of the matrix represent the variance of the corresponding parameters, and the non-diagonal elements represent the correlation between the corresponding two parameters. Figure 12 The anisotropy parameters calculated by inversion in Example 3 are compared with the exact values ​​in Table 6, where the error bars are Figure 11 The square roots of the diagonal elements of the covariance matrix in . Figure 13 is the average relative travel time error T in step 7 AV (H) Curve graph showing changes with preset H.

[0165] Example 4

[0166] In a specific embodiment of the present invention, Table 7 shows the elastic parameters of the monoclinic medium used in the implementation of the present invention. These elastic parameters comprise a 6×6 stiffness coefficient matrix, which is symmetric along the diagonal elements. Therefore, only the upper triangular matrix is ​​shown. In this embodiment, the reflector layer thickness is 0.7 km. The travel times of the horizontal reflection waves in the three-dimensional medium used in this embodiment were calculated using the ANRAY program package and used as the actual measured seismic data.

[0167]

[0168] Table 7 Model elastic parameters used in Example 4

[0169] Table 8 shows the reference P-wave velocity used in Example 4 of the present invention and the precise value of the anisotropy parameter calculated according to formula (20) based on the reference P-wave velocity and the elastic parameters in Table 7. The reference velocity can be defined as an arbitrary value. Changes in the reference velocity will only affect the anisotropy parameter itself and will not affect the accuracy of the longitudinal wave velocity described by the anisotropy parameter and the corresponding reflection process. Table 8 also shows the distribution range of the observation system in the present invention, where is the strike angle of the seismic line, and x is the length of the seismic line.

[0170]

[0171] Table 8 Observation system range and exact values ​​of anisotropy parameters in Example 4

[0172] Figure 14 The distribution range of the seismic survey line in Example 1 is drawn according to the data distribution range parameters in Table 8. Figure 15 is the covariance matrix calculated according to formula (33)-formula (35). The diagonal elements of the matrix represent the variance of the corresponding parameters, and the non-diagonal elements represent the correlation between the corresponding two parameters. Figure 16 The anisotropy parameters calculated by inversion in Example 4 are compared with the exact values ​​in Table 8, where the error bars are Figure 15 The square roots of the diagonal elements of the covariance matrix in . Figure 17 is the average relative travel time error T in step 7 AV (H) Curve graph showing changes with preset H.

[0173] Example 5

[0174] In a specific embodiment of the present invention, Table 9 shows the elastic parameters of the monoclinic medium used in the implementation of the present invention. These elastic parameters comprise a 6×6 stiffness coefficient matrix, which is symmetric along the diagonal elements. Therefore, only the upper triangular matrix is ​​shown. In this embodiment, the reflector thickness is 0.7 km. The travel time of the horizontal reflection wave in the three-dimensional medium used in this embodiment was calculated using the ANRAY program package and used as the actual measured seismic data.

[0175]

[0176]

[0177] Table 9 Model elastic parameters used in Example 5

[0178] Table 10 shows the reference P-wave velocity used in Example 5 of the present invention and the precise value of the anisotropy parameter calculated according to formula (20) based on the reference P-wave velocity and the elastic parameters in Table 9. The reference velocity can be defined as an arbitrary value. Changes in the reference velocity will only affect the anisotropy parameter itself and will not affect the accuracy of the longitudinal wave velocity described by the anisotropy parameter and the corresponding reflection process. Table 10 also shows the distribution range of the observation system in the present invention, where is the strike angle of the seismic line, and x is the length of the seismic line.

[0179]

[0180] Table 10 Observation system range and exact values ​​of anisotropy parameters in Example 5

[0181] Figure 18 is the distribution range of the seismic survey line in Example 1 drawn according to the data distribution range parameters in Table 10, Figure 19is the covariance matrix calculated according to formula (33)-formula (35). The diagonal elements of the matrix represent the variance of the corresponding parameters, and the non-diagonal elements represent the correlation between the corresponding two parameters. Figure 20 The anisotropy parameters calculated by inversion in Example 5 are compared with the exact values ​​in Table 10, where the error bars are Figure 19 The square roots of the diagonal elements of the covariance matrix in . Figure 21 This is a graph showing the average relative travel time error changing with the preset H.

[0182] Figure 22 This is a block diagram of an inversion device based on three-dimensional reflected P-wave anisotropy parameters provided in an embodiment of the present application. The device can be implemented as part or all of an electronic device through software, hardware, or a combination of both. Figure 22 As shown, the device includes:

[0183] A determination module is used to determine the spatial coordinate system of the seismic data, wherein the spatial coordinate system includes establishing a three-dimensional global coordinate system and a two-dimensional local coordinate system;

[0184] The local inversion module is used to invert the anisotropy parameters in a two-dimensional local coordinate system using the two-dimensional travel time data to obtain the first anisotropy parameter of the monoclinic medium;

[0185] A global inversion module is used to convert the first anisotropy parameter into a second anisotropy parameter in the three-dimensional global coordinate system based on a pre-established conversion relationship of anisotropy parameters between the three-dimensional global coordinate system and the two-dimensional local coordinate system, and to invert the global anisotropy parameter based on the second anisotropy parameter to obtain the second anisotropy parameter;

[0186] The calculation module is used to scan and calculate the formation thickness according to the second anisotropy parameter to obtain the target formation thickness of the underground medium.

[0187] The above description is only a preferred embodiment of the present application and is not intended to limit the scope of protection of the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application are included in the scope of protection of the present application.

[0188] The foregoing is merely a list of specific embodiments of the present application, intended to enable those skilled in the art to understand or implement the present application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application is not limited to the embodiments shown herein, but is intended to conform to the broadest scope consistent with the principles and novel features of the present application.

Claims

1. An inversion method based on anisotropic parameters of a monoclinic medium using three-dimensional P-wave reflection, characterized in that: include: Determining a spatial coordinate system for the seismic data, wherein the spatial coordinate system includes establishing a three-dimensional global coordinate system and a two-dimensional local coordinate system; Inverting the anisotropy parameter in the two-dimensional local coordinate system using the two-dimensional travel time data to obtain the first anisotropy parameter of the monoclinic medium; Based on a pre-established conversion relationship of anisotropy parameters between a three-dimensional global coordinate system and a two-dimensional local coordinate system, the first anisotropy parameter is converted into a second anisotropy parameter in the three-dimensional global coordinate system, and the global anisotropy parameter is inverted based on the second anisotropy parameter to obtain the second anisotropy parameter; A scanning calculation of the formation thickness is performed according to the second anisotropy parameter to obtain a calculation result.

2. The method according to claim 1, characterized in that The inversion of anisotropic parameters using the two-dimensional travel time data in the two-dimensional local coordinate system to obtain the first anisotropic parameter of the monoclinic medium includes: Construct an approximate formula for the variation of reflected P-wave travel time with offset distance in a two-dimensional local coordinate system; Constructing a first inversion determinant using the approximate formula for the variation of the reflected P-wave travel time with the offset; The first inversion determinant is solved using the least square method to obtain the first anisotropy parameter of the monoclinic medium.

3. The method according to claim 2, characterized in that The approximate formula for the change of reflected P-wave travel time with offset distance in the two-dimensional local coordinate system is as follows: Where x is the offset, T = T0(x) represents the travel time of the reflected P wave, and the anisotropic parameters include ε x '、ε z '、η y ',ε x '、ε z 'represents the difference between the P-wave phase velocity of the medium in the coordinate axis x1 and x3 directions and the reference velocity α; η y 'Represents the gradient of the P-wave velocity change in the medium in the coordinate plane (x1, x3).

4. The method according to claim 1, wherein The first inversion determinant is as follows: Gm=d; Where m = (ε x ',η y ',ε x ') T , G is an N×M coefficient matrix, N is the number of detectors, M is the number of A parameters involved, d is a vector, and the row vector in the matrix G is of the following form:

5. The method according to claim 1, wherein The first anisotropy parameter is converted to the second anisotropy parameter in the three-dimensional global coordinate system using the following formula: e' z =e z ; The second anisotropy parameter includes ε x ,ε y ,ε z ,η x ,η y ,η z , where ε x ,ε y ,ε z Respectively represent the difference between the P-wave phase velocity of the medium in the coordinate axes x1, x2, and x3 and the reference velocity α; η x , η y , η z represent the variation gradient of the P-wave velocity in the coordinate planes (x2, x3), (x1, x3), and (x1, x2), respectively. z represents the degree to which the P wave deviates from the axial direction in the plane (x1, x2); ξ 16 ,ξ 26 Represents the degree of deviation between the P wave direction and the fast and slow S wave directions in the plane (x1, x2).

6. The method according to claim 5, characterized in that The inversion of the global anisotropy parameter based on the second anisotropy parameter to obtain the second anisotropy parameter includes: The second inversion determinant is constructed based on the conversion relationship of anisotropic parameters between the 3D global coordinate system and the 2D local coordinate system; The second inversion determinant is solved using the least squares method to obtain a second anisotropy parameter.

7. The method according to claim 1, characterized in that The second inversion determinant is as follows: G1m1=d1; G2m2=d2; Among them, G1 is an N×6 coefficient matrix, N is the number of seismic lines distributed along different directions, G2 is an N×5 coefficient matrix, the number of rows N is the number of azimuths of the local coordinate system involved in the inversion, and the forms of vectors m1 and m2 are as follows: m1=(e x ,he y ,or z ,x z ,x 16 ,x 26 ) T ; m2=(η x ,or y ,or z ,x 16 ,x 26 ) T ; The forms of vectors d1 and d2 are as follows:

8. The method according to claim 1, characterized in that After performing inversion of a global anisotropy parameter based on the second anisotropy parameter to obtain the second anisotropy parameter, the method further includes: Calculate the residual between the model synthetic data and the actual seismic data during the inversion process; The inversion accuracy of the inversion of the global anisotropy parameters is determined by the residuals; When the reaction accuracy meets the preset conditions, a scanning calculation of the stratum thickness is performed based on the second anisotropy parameter and seismic data to obtain spatial information of the underground medium.

9. The method according to claim 1, characterized in that The scanning calculation of stratum thickness based on the second anisotropy parameter and seismic data to obtain spatial information of the underground medium includes: Get the thickness of the preset stratum; Calculating a first travel time of each detection point based on the second anisotropy parameter and the preset formation thickness; Obtaining the true second travel time of each detection point during the inversion process, and calculating the average relative travel time based on the first travel time and the second travel time; The stratum thickness with the smallest average relative travel time error is selected as the calculation result.

10. An inversion device based on three-dimensional reflected P-wave anisotropy parameters, characterized in that: include: A determination module, configured to determine a spatial coordinate system of seismic data, wherein the spatial coordinate system includes establishing a three-dimensional global coordinate system and a two-dimensional local coordinate system; A local inversion module is used to invert anisotropic parameters in the two-dimensional local coordinate system using the two-dimensional travel time data to obtain a first anisotropic parameter of the monoclinic medium; a global inversion module, configured to convert the first anisotropy parameter into a second anisotropy parameter in the three-dimensional global coordinate system based on a pre-established conversion relationship of anisotropy parameters between the three-dimensional global coordinate system and the two-dimensional local coordinate system, and perform an inversion of the global anisotropy parameter based on the second anisotropy parameter to obtain the second anisotropy parameter; The calculation module is used to scan and calculate the formation thickness according to the second anisotropy parameter to obtain the target formation thickness of the underground medium.

Images