A joint focusing inversion method for heavy earthquakes based on unstructured grids

Abstract

The present invention relates to the field of joint inversion, and specifically to a gravity-seismic joint focusing inversion method under unstructured grid division, comprising the following steps: performing unstructured tetrahedral grid division on the underground space; performing forward modeling on gravity data and seismic data respectively to obtain forward gravity data and forward seismic data; establishing an objective function for the forward gravity data and forward seismic data using a constraint function of minimum support of a fractal index; and obtaining density joint inversion results and velocity joint inversion results. The present invention can obtain higher-resolution joint inversion results in practical applications. Compared with traditional focused joint inversion ideas, the present invention can more finely characterize the distribution range of physical properties and complex stratum boundaries, so that gravity-seismic joint inversion under unstructured tetrahedral grid division with higher inversion resolution and field source characterization capabilities has better development in fields such as fine structure exploration.

Description

Technical Field

[0001] The present invention relates to the field of joint inversion, in particular to a heavy earthquake joint focusing inversion method under unstructured grid subdivision. Background Art

[0002] Joint inversion of heavy earthquake data is a method of inverting the distribution of underground physical properties using observation data. It requires dividing the underground space to be inverted into closely arranged grid cells, and establishing a kernel matrix through the forward formula of the grid cell at the surface observation point under unit physical properties. Then, a linear equation system is established using the measured data at the observation point, the kernel matrix, and the unknown physical properties at each grid cell. The physical property values ​​at the underground grid cell are obtained by solving the linear equation system.

[0003] In joint inversion of geophysical data, two common meshing methods are structured hexahedral meshing and unstructured tetrahedral meshing. Unstructured tetrahedral meshing can better depict the undulating terrain and underground strata distribution, and has higher inversion resolution. Therefore, to achieve better inversion results, unstructured tetrahedral meshing is now preferred in joint inversion.

[0004] Currently, the most common joint inversion technique is cross-gradient joint inversion, while focused inversion is a method that can obtain higher-resolution inversion results. However, this method is more commonly used in separate geophysical inversions using structured hexahedral grids. For joint inversion, a joint minimum focusing term has been added to the joint inversion of electromagnetic and gravity gradient data, and a joint minimum entropy constraint has been adopted for gravity-magnetic joint inversion. These methods are still only applied to structured hexahedral grids and suffer from insufficient resolution of the inversion results and dependence on multiple parameters. Therefore, existing focusing strategies cannot be used in focused joint inversion using unstructured tetrahedral grids. This severely limits the application of joint inversion methods using gravity-seismic data using unstructured tetrahedral grids in resource exploration and the ability to accurately characterize underground structures. Summary of the Invention

[0005] One approach to geophysical focused joint inversion on existing structured hexahedral grids 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 goal of focusing the result. However, the focusing ability of the joint minimum support function is insufficient, and higher-resolution inversion results cannot be obtained. Another approach 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, because it introduces too many parameters that affect the focusing results as variables, it is difficult to determine the problem in solving different practical problems and is not very practical. Currently, focused joint inversion is still only used in single geophysical inversion on structured hexahedral grids and a few joint inversions. Focused joint inversion methods on unstructured tetrahedral grids are lacking and computationally complex. Therefore, there is a greater demand for how to obtain higher-resolution joint inversion results on unstructured tetrahedral grids. This invention proposes for the first time a constraint function based on the minimum support of the fractal index, which is effectively applied to the joint focus inversion of gravity and seismic data on unstructured tetrahedral grids, improving the resolution of the joint inversion results. It is mainly used in resource exploration, underground physical property distribution, and detailed characterization of stratigraphic structure. The invention provides the following technical solutions:

[0006] A heavy-seismic joint focusing inversion method based on unstructured grid subdivision includes the following steps:

[0007] The first step is to use unstructured tetrahedral mesh to divide the underground space;

[0008] The second step is to forward model the gravity data and seismic data to obtain forward gravity data and forward seismic data;

[0009] 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:

[0010] ;

[0011] ;

[0012] ,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 ability of different types of focus constraints, Representative physical properties, is the focusing factor, usually selected , It 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 it at the point where the slope of the curve is the largest. The optimal focusing factor ;

[0013] The fourth step is to obtain the density joint inversion results and the velocity joint inversion results.

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

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

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

[0017] As a further solution of the present invention: the expression of the type function is:

[0018] ;

[0019] ;

[0020] ;

[0021] ;

[0022] ; represents the interpolation coefficient, represents the complex conjugate, represents the fundamental frequency, Is a random number and guarantees , the center point coordinates are , non-zero variables is the scaling factor, It is a function of different orders. The calculation formula of the physical property gradient is:

[0023] ,

[0024] ,

[0025] , ,

[0026] .

[0027] Compared with the prior art, the present invention has the following beneficial effects:

[0028] In practical applications, the present invention can obtain higher-resolution joint inversion results. Compared with the traditional focused joint inversion approach, the present invention can more precisely characterize the distribution range of physical properties and complex stratigraphic boundaries, enabling the joint inversion of gravity seismic data under unstructured tetrahedral meshing with higher inversion resolution and field source characterization capabilities to have better development in fields such as fine structural exploration. The present invention is the first to implement joint focused inversion of gravity seismic data under unstructured tetrahedral meshing. First, we give a threshold and determine the ratio of the number of nodes greater than the threshold to the actual volume of the model as the quantitative resolution. Through simulation experiments, it is verified that the resolution of the joint inversion results of gravity seismic data under unstructured tetrahedral meshing is improved by 7%. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 Flowchart of the heavy-seismic joint focusing inversion method under unstructured grid subdivision in an embodiment of the present invention.

[0030] Figure 2 Comparison of the accuracy of different order functions in the combined focusing inversion method for heavy seismicity on an unstructured grid in an embodiment of the present invention. (a) shows the density model of the theoretical cosine function; (b) shows the velocity model of the theoretical cosine function; (c) shows the box plots of different order functions; (d) shows the error between the calculated result of the first-order function and the true model; (e) shows the error between the calculated result of the second-order function and the true model; and (f) shows the error between the calculated result of the third-order function and the true model.

[0031] Figure 3 3 is a comparison chart of the theoretical focusing capabilities of the gravity-seismic joint focusing inversion method under unstructured grid subdivision in an embodiment of the present invention and the traditional method.

[0032] Figure 4Comparison of results from different joint focusing inversion methods for tilted prisms. (a) shows the forward gravity modeling results and model distribution for the tilted prism model; (b) shows the wavefield simulation results; (c) shows the density results from the joint focusing inversion of gravity seismic events using the minimum support algorithm on an unstructured grid; (d) shows the velocity results from the joint focusing inversion of gravity seismic events using the minimum support algorithm on an unstructured grid; (e) shows the density results from the joint focusing inversion of gravity seismic events using the algorithm of the present invention on an unstructured grid; and (f) shows the velocity results from the joint focusing inversion of gravity seismic events using the algorithm of the present invention on an unstructured grid.

[0033] Figure 5 Comparison of results from different joint focusing inversion methods for a tilted formation model. (a) shows the forward gravity modeling results and model distribution of the formation model; (b) shows the wavefield simulation results; (c) shows the density results from the combined minimum support algorithm for gravity-seismic joint focusing inversion on an unstructured grid; (d) shows the velocity results from the combined minimum support algorithm for gravity-seismic joint focusing inversion on an unstructured grid; (e) shows the density results from the combined minimum support algorithm for gravity-seismic joint focusing inversion on an unstructured grid; and (f) shows the velocity results from the combined minimum support algorithm for gravity-seismic joint focusing inversion on an unstructured grid. DETAILED DESCRIPTION

[0034] The following will provide a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0035] The objective function formula of the joint minimum support focusing inversion method under the traditional structure hexahedron subdivision is:

[0036] ;

[0037] ;

[0038] ;

[0039] ;

[0040] in are the data fitting term and the regularization term. and These are two types of physical property information corresponding to each subdivision unit. N is the number of measured data points, and is the observation data, and is the data for forward calculation, Yes and The error is weighted to account for the data differences. W m is the depth weight function, which is used to balance the weights of different depth subdivision units. m 1ref and m 2ref It is a reference model. is the joint minimum support focus term, where Yes, the initial model. is a minimum value to prevent singularities.

[0041] This paper proposes a joint focusing inversion technique for unstructured tetrahedral grid re-seismic with minimum support of fractal index, constructs a new fractal index minimum support focusing constraint term and adds it to the objective function. The specific form is as follows:

[0042] ;

[0043] ;

[0044] ;

[0045] ;

[0046] in are the data fitting term and the regularization term. and is the density and velocity corresponding to each subdivision unit, 、 The present invention introduces the classification weight function 、 , which has the ability to balance different types of focusing constraints, is firstly The function quickly obtains the range of the target body. As the number of iterations increases, the The weight of the structure is used to enhance the structural consistency to define the boundary of the target body. In order to prevent excessive structural constraints from causing distortion of the inversion results or incorrect data fitting, the present invention reduces the influence of structural constraints, enhances data consistency, and ultimately achieves the purpose of data fitting. is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation, usually chosen .

[0047] The calculation of physical property gradients under unstructured tetrahedral grids is included. However, the traditional difference and FFT methods cannot be used for calculation under tetrahedral grids, and the calculation is complex and time-consuming, which limits its application. Research is carried out on the calculation of physical property gradients under unstructured tetrahedral grids, and a NUFFT fast calculation method based on second-order finite difference functions is proposed.

[0048] ;

[0049] is the data that meets the FFT operation conditions, M is the number of sampling points, is the inversion result, It is a type function, and its specific expression is as follows:

[0050] ;

[0051] ;

[0052] ;

[0053] ;

[0054] ;

[0055] represents the interpolation coefficient, represents the complex conjugate, represents the fundamental frequency, Is a random number and guarantees The center point coordinates are , non-zero variables is the scaling factor, Different order functions require different number of points to participate in the calculation under unstructured tetrahedral mesh. The first order function is determined by four coefficients, so it needs 4 points around to participate in the calculation. Similarly, the second order function requires 10 points, and the third order function requires 20 points. The higher the order of the function, the higher the accuracy of the data, but the longer the calculation time under unstructured tetrahedral mesh. Therefore, it is meaningful to select a few order functions under the premise of ensuring accuracy. The calculation formula for the physical property gradient is as follows:

[0056] ,

[0057] ,

[0058] , ,

[0059] .

[0060] The specific implementation of the present invention is described in detail below with reference to specific embodiments.

[0061] See Figure 1-Figure 3 , using the box plot method to compare the capabilities of different order functions in terms of accuracy and computational efficiency, the effect is compared. Figure 2 shown. Figure 2 a, 2b are the density model and velocity model of the theoretical cosine function. Figure 2 d, 2e, and 2f are the errors between the results of the first-order, second-order, and third-order functions and the true model. It can be seen that the first-order accuracy is low, while the second-order and third-order accuracy are good. At the same time, the box plot is compared. Figure 2 The results of c show that the accuracy of the second-order and third-order functions is close, but due to the different amount of calculation required to calculate a single point, the time for the first-order calculation of a single point is 0.11 seconds, the second-order is 0.25 seconds, and the third-order is 0.43 seconds. In addition, as the number of underground grids increases, the amount of calculation increases. Finally, considering the calculation accuracy and efficiency, it is more reasonable to determine the second-order function. Theoretical focusing ability comparison, effect comparison Figure 3 As shown, Figure 3 The horizontal axis is the physical property of the real model m Physical properties of reference models m 0 difference (when there is no reference model, m 0=0), the vertical axis is the value of different focusing items. When the horizontal axis is the same, the method proposed in the present invention has a larger value than the traditional joint minimum support method. Therefore, the method proposed in the present invention has better focusing ability than the traditional joint minimum support method.

[0062] The method of the present invention is used to calculate the combined inversion of two underground inclined prisms and the inclined stratum model. The inversion results are as follows: Figure 4 and Figure 5 As 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 traditional methods, Figure 4 e, 4f are the density and velocity distributions obtained by the method of the present invention, Figure 4 It can be clearly seen 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 traditional methods, Figure 5 5e, 5f are the density and velocity distributions obtained by the method of the present invention, Figure 5It can be clearly seen that the present invention can obtain higher resolution results, the distribution of the results and the physical property recovery values ​​are closer to the real ones, and 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 cutoff parameters and count the proportion of the number greater than the physical property cutoff within the target range. In the inclined prism model, we set the density to be greater than 0.2 g / cm 3 The velocity is greater than 1950 m / s as the cutoff parameter. 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 to be greater than 0.5 g / cm 3 The speed is greater than 2600 m / s as the cutoff parameter. Compared with the traditional method, the resolution of the joint inversion result of the present invention is improved by 7.4%.

[0063] In addition, it should be understood that although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A joint focusing inversion method for heavy earthquakes based on unstructured grid subdivision, characterized by: The following steps are involved: 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 gravity data and forward 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, There are physical property gradients in unstructured tetrahedral grids. is a classification weight function that balances the ability of different types of focus constraints, Representative physical properties, is the focusing coefficient, is the focusing factor; The fourth step is to obtain the density joint inversion results and the velocity joint inversion results.

2. The method for combined focusing inversion of heavy-seismic data under unstructured grid subdivision according to claim 1 is characterized in that: described , is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation.

3. The heavy-seismic joint focusing inversion method based on unstructured grid subdivision according to claim 1 or 2 is characterized in that: described , is the regularization parameter, represents the number of iterations in the inversion process, is the standard deviation.

4. The method for combined focusing inversion of heavy-seismic data under unstructured grid subdivision according to claim 1 is characterized in that: described The physical property gradient existing in the unstructured tetrahedral grid is calculated using the NUFFT fast calculation method of the second-order finite difference 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, It is a type function.

Images