Gravity-earthquake combined focusing inversion method under unstructured grid subdivision

Abstract

The invention relates to the field of joint inversion, in particular to a gravity-earthquake joint focusing inversion method under unstructured mesh generation, which comprises the following steps of: performing unstructured tetrahedral mesh generation on an underground space by using an unstructured tetrahedral mesh; forward modeling is carried out on the gravity data and the seismic data to obtain forward modeling gravity data and forward modeling seismic data; establishing an objective function for the forward modeling gravity data and the forward modeling seismic data by using a constraint function supported by the minimum parting index; and obtaining a density joint inversion result and a speed joint inversion result. According to the method, a joint inversion result with higher resolution can be obtained in practical application, and compared with a traditional focusing joint inversion thought, the physical property distribution range and the complex stratum boundary line can be depicted more finely; therefore, the heavy-seismic joint inversion under the non-structural tetrahedral mesh generation with higher inversion resolution and field source depicting capability has better development in the fields of fine structure exploration and the like.

Description

Technical Field

[0001] The present invention relates to the field of joint inversion, and specifically to a gravity-seismic joint focusing inversion method under unstructured grid meshing. Background Art

[0002] Gravity-seismic data joint inversion is a method for inversely determining the subsurface physical property distribution using observed data. It is necessary to divide the subsurface space to be inverted into closely arranged grid cells, establish a kernel matrix through the forward modeling formula of grid cells with unit physical properties at surface observation points, and then establish a linear equation system using the measured data at observation points, the kernel matrix, and the physical property unknowns at each grid cell. The physical property values at subsurface grid cells are obtained by solving the linear equation system.

[0003] In geophysical data joint inversion, there are two commonly used grid meshing methods, namely structured hexahedral grids and unstructured tetrahedral grids. Using unstructured tetrahedral grid meshing can better depict the undulating terrain and subsurface strata distribution, and has higher inversion resolution. Therefore, for better inversion results, current joint inversion tends to use unstructured tetrahedral grid meshing.

[0004] Currently, the general joint inversion technical means is cross-gradient joint inversion, and the focusing inversion method is a method that can obtain a higher-resolution inversion result. However, this method is more applied to the individual inversion of structured hexahedral grid geophysics. For joint inversion, currently, a joint minimum focusing term is added to the joint inversion of electromagnetic and gravity gradient data, and a joint minimum entropy constraint method is adopted for gravity-magnetic joint inversion. These methods are still only applied to structured hexahedral grid meshing, and there are problems such as insufficient resolution of inversion results and dependence on multiple parameters. Therefore, the existing focusing strategies cannot be used in the focusing joint inversion under unstructured tetrahedral grid meshing, which severely limits the application of the gravity-seismic data joint inversion method under unstructured tetrahedral grid meshing in resource exploration and the ability to finely depict the subsurface structure. Summary of the Invention

[0005] One way of focusing joint inversion of geophysics under existing structured hexahedral grid division is to introduce a single joint minimum support function, which is proportional to the model volume. Minimizing the objective function is equivalent to continuously reducing the volume of the inversion result, thereby achieving the purpose of focusing the result. However, the focusing ability of the joint minimum support function is insufficient, and it is impossible to obtain higher resolution inversion results. Another way is to add joint minimum entropy. This method enhances the similarity between parameters of different types of models, thereby obtaining high-resolution inversion results. However, since it introduces too many parameters that affect the focusing results as variables, there are difficult problems that are difficult to determine in solving different practical problems, and it does not have good practicality. At present, focused joint inversion is still only used in single inversion of structured hexahedral grid geophysics and a few joint inversions. The focused joint inversion method under unstructured tetrahedral grid is lacking and the calculation is complicated. Therefore, there is a greater demand for how to obtain higher resolution joint inversion results under unstructured tetrahedral grid. The present invention proposes for the first time a constraint function based on the minimum support of the fractal index, which is effectively applied to the gravity-seismic joint focusing inversion of unstructured tetrahedral grids, improving the resolution of the gravity-seismic joint inversion results, and is mainly used in resource exploration, underground physical property distribution and fine characterization of stratum structure. The present invention provides the following technical solutions: A heavy-seismic joint focusing inversion method based on unstructured grid subdivision includes the following steps: The first step is to use unstructured tetrahedral mesh to divide the underground space; The second step is to forward model the gravity data and seismic data to obtain forward modeled gravity data and forward modeled seismic data; In the third step, the objective function is established by using the constraint function of the minimum support of the classification index for the forward gravity data and the forward seismic data. The objective function is in the following form: ; ; ,in are the data fitting term and the regularization term, and is the density and velocity corresponding to each subdivision unit, are different types of focus constraints constructed, is a classification weight function that balances the capabilities of different types of focusing constraints, Representative physical properties, is the focusing factor, usually chosen , is the focusing factor. The principle of selecting the focusing factor for different focusing items is to use a method similar to the L curve and select at the point where the slope of the curve is the largest. The optimal focusing factor ; Step 4: Obtain the density joint inversion result and the velocity joint inversion result.

[0006] As a further solution of the present invention: , is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation, usually selected .

[0007] As a further solution of the present invention: , is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation, usually selected .

[0008] As a further solution of the present invention: Calculate the physical property gradient under the unstructured tetrahedral mesh in using the NUFFT fast calculation method of the second-order finite-difference type function. The formula is as follows: is the data that meets the FFT operation conditions, M is the number of sampling points, is the inversion result, is the type function.

[0009] As a further solution of the present invention: The expression of the type function is: ; ; ; ; ; represents the interpolation coefficient, represents the complex conjugate, represents the fundamental frequency, is a random number and ensures , the center point coordinates are , the non-zero variable is the scaling factor, is the function of different orders. At this time, the calculation formula of the physical property gradient is: , , , , .

[0010] Compared with the prior art, the beneficial effects of the present invention are: In practical applications, the present invention can obtain a joint inversion result with higher resolution. Compared with the traditional focused joint inversion idea, the present invention can more precisely depict the physical property distribution range and complex formation boundaries, enabling the gravity-seismic joint inversion under unstructured tetrahedral mesh division with higher inversion resolution and field source characterization ability to have better development in fields such as fine structure exploration. The present invention realizes the gravity-seismic joint focused inversion for the first time under unstructured tetrahedral mesh division. First, we set a threshold and determine the ratio of the number of nodes greater than the threshold to the true volume of the model as the quantitative resolution. Through simulation experiments, it is verified that the resolution of the gravity-seismic joint inversion result under unstructured tetrahedral mesh division is increased by 7%. Brief Description of the Drawings

[0011] Figure 1 It is a flow chart of the gravity-seismic joint focused inversion method under unstructured mesh division in an embodiment of the present invention.

[0012] Figure 2 It is a comparison chart of the accuracies of different order type functions in the gravity-seismic joint focused inversion method under unstructured mesh division in an embodiment of the present invention. Among them, (a) is the density model of the theoretical cosine function; (b) is the velocity model of the theoretical cosine function; (c) is the box plot of different order type functions; (d) is the error between the result calculated by the first-order type function and the true model; (e) is the error between the result calculated by the second-order type function and the true model; (f) is the error between the result calculated by the third-order type function and the true model.

[0013] Figure 3 It is a comparison chart of the theoretical focusing capabilities between the gravity-seismic joint focused inversion method under unstructured mesh division in an embodiment of the present invention and the traditional method.

[0014] Figure 4 It is a comparison chart of the results of different joint focused inversion methods for an inclined prism. Among them, (a) is the gravity forward inversion result and model distribution of the inclined prism model; (b) is the wave field simulation result; (c) is the density result of the gravity-seismic joint focused inversion under unstructured mesh division using the joint minimum support algorithm; (d) is the velocity result of the gravity-seismic joint focused inversion under unstructured mesh division using the joint minimum support algorithm; (e) is the density result of the gravity-seismic joint focused inversion under unstructured mesh division using the algorithm of the present invention; (f) is the velocity result of the gravity-seismic joint focused inversion under unstructured mesh division using the algorithm of the present invention.

[0015] Figure 5Comparison chart of results of different joint focusing inversion methods for an inclined formation model. Among them, (a) is the gravity forward modeling result of the formation model and the model distribution; (b) is the wave field simulation result; (c) is the density result of joint gravity-seismic focusing inversion with unstructured grid meshing using the joint minimum support algorithm; (d) is the velocity result of joint gravity-seismic focusing inversion with unstructured grid meshing using the joint minimum support algorithm; (e) is the density result of joint gravity-seismic focusing inversion with unstructured grid meshing using the algorithm of the present invention; (f) is the velocity result of joint gravity-seismic focusing inversion with unstructured grid meshing using the algorithm of the present invention. Detailed implementation manners

[0016] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without making creative efforts shall fall within the protection scope of the present invention.

[0017] The objective function formula of the joint minimum support focusing inversion method under traditional structured hexahedron meshing is: ; ; ; ; Among them is the data fitting term and the regularization term. and are two kinds of physical property information corresponding to each meshing unit, N is the number of measured data points, and are the observed data, and are the forward calculation data, is and is the error, which weights the data difference. W m is the depth weight function, which is used to balance the weights of meshing units at different depths, m 1ref and m 2ref are the reference models. is the joint minimum support focusing term, where is the initial model of, is a minimum value to prevent the occurrence of singularities.

[0018] The present invention proposes a non - structured tetrahedral mesh joint seismic and focusing inversion technology supported by the minimum fractal index. A new focusing constraint term supported by the minimum fractal index is constructed and added to the objective function, and its specific form is as follows: ; ; ; ; where is the data fitting term and the regularization term. and are the density and velocity corresponding to each subdivision unit. , are different types of focusing constraint terms constructed. The present invention introduces the fractal weight functions , , which have the ability to balance different types of focusing constraint terms. During the inversion process, first, the range of the target body is quickly obtained through the function. As the number of iterations increases, the weight of is increased to enhance the structural consistency to delimit the boundary of the target body. To prevent the inversion result from being distorted or incorrect data fitting caused by excessive structural constraints, the present invention weakens the influence of structural constraints and enhances data consistency to ultimately achieve the purpose of data fitting. is the regularization parameter, represents the number of iterations during the inversion process, is the standard deviation, usually selected as .

[0019] contains the calculation of the physical property gradient under the non - structured tetrahedral mesh. Since the traditional difference and FFT methods cannot be used for operations under the tetrahedral mesh, and due to problems such as its complex calculation and time - consuming calculation, its application is limited. Therefore, the present invention conducts research on the calculation of the physical property gradient under the non - structured tetrahedral mesh in and proposes a fast NUFFT calculation method based on the second - order finite - difference type function.

[0020] ; is the data that meets the FFT operation conditions, M is the number of sampling points, is the inversion result, is the type function, and its specific expression form is as follows: ; ; ; ; ; represents the interpolation coefficient, denotes the complex conjugate, represents the fundamental frequency, is a random number and ensures that . The center point coordinates are , and the non-zero variable is the scaling factor, is a function of different orders. For non-structured tetrahedral meshes, the number of points required for different order shape functions to participate in the calculation is different. The first-order shape function is determined by four coefficients, so 4 surrounding points are required to participate in the calculation. Similarly, the second-order shape function requires 10 points, and the third-order shape function requires 20 points. The higher the order of the shape function, the higher the data accuracy, but the longer the calculation time for non-structured tetrahedral meshes. Therefore, it is of research significance to select the order of the shape function on the premise of ensuring accuracy. The calculation formula for the physical property gradient is as follows: , , , , .

[0021] The following describes the specific implementation of the present invention in detail with reference to specific embodiments.

[0022] Refer to Figures 1 - 3 , and use the box plot method to compare the capabilities of different order shape functions in terms of accuracy and computational efficiency. The comparison results are as Figure 2 shown. Figure 2 a, 2b are the density model and velocity model of the constructed theoretical cosine function, Figure 2 d, 2e, 2f are the errors between the calculation results of the first-order, second-order, and third-order shape functions and the true model. It can be seen that the first-order accuracy is relatively low, and the second-order and third-order accuracies are better. At the same time, comparing the results of the box plot Figure 2 c, it can be seen that the second-order and third-order shape function accuracies are close, but due to the different amounts of computation required for calculating a single point, the time for the first-order to calculate a single point is 0.11 seconds, the second-order is 0.25 seconds, and the third-order is 0.43 seconds. And as the underground dissection grid becomes more, the amount of computation becomes larger. Finally, considering both computational accuracy and computational efficiency, it is determined that the second-order shape function is more reasonable. Conduct a comparison of the theoretical focusing ability. The comparison results are as Figure 3 shown, Figure 3 in which the abscissa is the difference between the physical property of the true model m and the physical property of the reference model m 0 (when there is no reference model, m0 = 0), the ordinate is the value of different focusing terms. At the same abscissa, the method proposed by the present invention has a larger value than the traditional combined minimum support method. Therefore, the method proposed by the present invention has better focusing ability compared with the traditional combined minimum support method.

[0023] By using the method of the present invention, the joint gravity-seismic inversion of two underground inclined prisms and an inclined formation model is calculated respectively. The inversion results are as Figure 4 and Figure 5 shown. Figure 4 a is the inclined prism and its corresponding gravity anomaly data, Figure 4 b is the seismic wave field record data, Figure 4 c, 4d are the density and velocity distributions obtained by the traditional method, Figure 4 e, 4f are the density and velocity distributions obtained by the method of the present invention. It can be clearly seen from Figure 4 that the present invention can obtain results with higher resolution, and the distribution of the results and the physical property recovery values are closer to the real ones. Figure 5 a is the formation model and its corresponding gravity anomaly data, Figure 5 b is the seismic wave field record data, Figure 5 c, 5d are the density and velocity distributions obtained by the traditional method, Figure 5 e, 5f are the density and velocity distributions obtained by the method of the present invention. It can be clearly seen from Figure 5 that the present invention can obtain results with higher resolution, and the distribution of the results and the physical property recovery values are closer to the real ones. The tests of different models can better verify the universality of the method of the present invention. In summary, the method of the present invention can obtain the position and boundary of the target body more accurately, and the physical property recovery is closer to the real value. In order to more quantitatively compare the improvement in resolution between the present invention and the traditional method, we first set the physical property truncation parameter and count the proportion of the number greater than the physical property truncation within the target range. In the inclined prism model, we set the density greater than 0.2 g / cm 3 and the velocity greater than 1950 m / s as the truncation parameters. Compared with the traditional method, the resolution of the joint inversion result of the present invention is improved by 7.2%; in the formation model, we set the density greater than 0.5 g / cm 3 and the velocity greater than 2600 m / s as the truncation parameters. Compared with the traditional method, the resolution of the joint inversion result of the present invention is improved by 7.4%.

[0024] In addition, it should be understood that although this specification is described according to the embodiments, not every embodiment only contains an independent technical solution. This narrative way of the specification is only for clarity. Those skilled in the art should regard the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A combined gravity and seismic joint focusing inversion method under unstructured grid meshing, characterized in that Including the following steps: First step, perform unstructured tetrahedral mesh generation on the underground space using unstructured tetrahedral meshes; Second step, respectively perform forward modeling on the gravity data and seismic data to obtain the forward modeled gravity data and forward modeled seismic data; Third step, establish an objective function using the constraint function supported by the minimum fractal index for the forward modeled gravity data and forward modeled seismic data. The form of the objective function is as follows: ； ； , where is the data fitting term and the regularization term, and are the density and velocity corresponding to each sub - division unit, is the constructed focusing constraint term of different types, There is a physical property gradient under the unstructured tetrahedral mesh, is the classification weight function for balancing the ability of different types of focusing constraint terms, represents the physical property type, is the focusing coefficient, is the focusing factor; Fourth step, obtain the density joint inversion result and the velocity joint inversion result.

2. The combined gravity and seismic joint focusing inversion method under unstructured grid meshing according to claim 1, wherein The , is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation.

3. The combined gravity and seismic joint focusing inversion method under unstructured grid meshing according to claim 1 or 2, characterized in that, The said , is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation.

4. The combined gravity and seismic joint focusing inversion method under unstructured grid meshing according to claim 1, characterized in that The The physical property gradient existing under the unstructured tetrahedral mesh is calculated by using the NUFFT fast calculation method of the second-order finite-difference type function, and the formula is as follows: is the data that meets the FFT operation conditions, M is the number of sampling points, is the inversion result, is the type function.

5. The combined gravity and seismic joint focusing inversion method under unstructured grid meshing according to claim 4, wherein The expression of the said type function is: ； ； ； ； ； represents the interpolation coefficient, denotes the complex conjugate, represents the fundamental frequency, is a random number and ensures that , the center point coordinates are , non - zero variable is the scaling factor, are functions of different orders.

Images