Model feature driven adaptive sparse constraint inversion method for gravity data

Abstract

The present application is suitable for the field of mineral geophysical exploration and application, and provides a gravity data adaptive sparse constraint inversion method based on model feature driving, which comprises: grid processing on collected gravity field data, generating regular grid data and completing three-dimensional discrete grid sectioning of underground space; constructing a three-dimensional gravity sparse constraint inversion objective function; model local feature identification, calculating model smoothness difference; model feature partition constraint, dividing model area according to smoothness difference; adaptive order dynamic optimization; iterative inversion solution, repeatedly executing the above order dynamic optimization until the convergence condition is met, and outputting the final three-dimensional density model. The present application realizes the sparse constraint inversion method of gravity data by combining model feature driven partition and iterative optimization, improves the inversion resolution and accuracy of the geological model, and enhances the applicability of the method technology to the interpretation of complex geological structures.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of mineral geophysical exploration and application, and particularly to a gravity data adaptive sparse constraint inversion method based on model feature driving. BACKGROUND

[0002] Gravity inversion is a mathematically ill-posed problem of observing gravity field data to deduce the underground density structure, which has the characteristics of non-unique solution and sensitivity to noise, etc. Therefore, a regularization method is usually introduced to constrain the target equation to obtain a reasonable density parameter model solution. The objective function of the regularization inversion of gravity data is usually composed of two parts: a data fitting term and a model constraint term. The data fitting term usually uses the 2-norm to measure the difference between the gravity observation data and the density model prediction data, and the model constraint term is to impose prior information constraints on the density parameters to obtain the optimal solution in the mathematical and physical sense. L

[0003] The diversified occurrence of underground geological structures will produce characteristics such as sudden change or continuous change of field source density parameters, so the field source structure can be simply classified into sharp edge structure and smooth structure. Sparse constraint inversion is to constrain the model solution space by using the L p -norm form, so as to represent the sharp edge and smooth structure characteristics of the field source. Mathematically, L p -norm constraint is defined as the sum of the absolute values of the elements of the model parameter vector to the power of p . On the one hand, the existing gravity data sparse constraint inversion based on the L p -norm only selects a fixed order p value for the whole model space, which cannot flexibly cope with the complex scenarios of coexistence of rock ore bodies with sharp edge characteristics, intrusive rock contact zones, and continuous sedimentary strata and altered rocks with smooth structure characteristics. p If the value is too small (such as p ≤1), the continuity of the smooth area of the model will be destroyed, p if the value is too large (such as p ≥1), the boundary of the block area with sharp edge characteristics will be blurred. On the other hand, the selection of the order p in sparse constraint inversion depends heavily on reliable prior information and trial and error, and the optimal value is often optimized through subjective experience or repeated experiments, which is difficult to quickly apply to inversion interpretation of complex geological structures. Systematic adaptive mechanism and stable and reliable sparse constraint inversion calculation scheme limit the resolution and accuracy of gravity data inversion in fine description of complex geological structures. SUMMARY

[0004] ​In view of the deficiencies of the prior art, the purpose of the present application is to propose a model feature driven gravity data adaptive sparse constraint inversion method to solve the technical problems mentioned in the background art.

[0005] The model feature driven gravity data adaptive sparse constraint inversion method comprises the following steps:

[0006] S1, grid processing is performed on the collected gravity field data to generate regular grid data and complete three-dimensional discrete grid dissection of the underground space;

[0007] S2, a three-dimensional gravity sparse constraint inversion objective function is constructed;

[0008] S3, model local feature recognition is performed, and the smoothness difference of the model is calculated;

[0009] S4, model feature partition constraint is performed, and the model area is divided according to the smoothness difference;

[0010] S5, adaptive order dynamic optimization is performed;

[0011] S6, iterative inversion solving is performed, the order dynamic optimization of S5 is repeatedly executed until the convergence condition is met, and the final three-dimensional density model is output.

[0012] Preferably, in step S2, the three-dimensional gravity sparse constraint inversion objective function has the following form:

[0013] Formula 1: ;

[0014] wherein, is a gravity field data fitting term for measuring the fitting degree of the observed gravity field data and the gravity field data predicted by the obtained density model, is a model constraint term for applying prior information constraint conditions to the density model to obtain a reasonable solution of the model space, represents the square operation of the vector l 2 norm, represents the lp norm of the vector p power operation, μ is a regularization parameter, p is a norm order, W d is a data weight matrix, which is usually composed of the inverse of data error, W m is a model weighting function, G is a forward kernel matrix of the gravity field data, d represents gravity field observation data, m is a model parameter representing the density distribution of the underground medium respectively;

[0015] To solve the zero-point non-differentiability, an approximate form of the norm is used lp instead of the norm lp which is strictly defined, to construct the model constraint term φ m An approximate form of (m) is as follows:

[0016] Equation 2: ;

[0017] wherein, m i represents each discrete grid cell of the subsurface space density model, ε is a small threshold value, ranging from 10 -1 ~10 -6 , to avoid the denominator equal to zero, and to enhance numerical stability.

[0018] Preferably, in step S3, the model local feature recognition, the content of calculating the model smoothness difference includes:

[0019] First, by fixing the order p =1, two sets of threshold values of different orders of magnitude difference are given ε , respectively, to complete the gravity field data sparse constraint inversion calculation, and obtain two independent three-dimensional density solution models m 1 and m 2, since the change rate of the model smoothness is strongly dependent on the inherent characteristics of the corresponding area in the true model, the blocky area that retains sharp edge features shows higher sensitivity to the threshold value ε ; therefore, by defining the smoothness of the two inversion model spaces L i , the differentiation of different characteristic areas is realized by the smoothness difference Δ L , and the smoothness measurement operator L i is defined as the sum of the second-order derivatives of the model space, and Δ L is the smoothness difference of the inversion model under different threshold values ε , which is as follows:

[0020] Equation 3: ;

[0021] wherein, m i represents each discrete grid cell of the subsurface space density model, i represents two sets of threshold values of different orders of magnitude difference ε independent inversion corresponding to the number, x , y , zrespectively represent three orthogonal directions in the Cartesian coordinate system.

[0022] Preferably, in step S4, the model feature partition constraint, which divides the content of the model region according to the smoothness difference, includes:

[0023] According to the size of the smoothness difference Δ L , the blockiness and smoothness of the model space are indicated, and the region with a larger smoothness difference (greater than the statistical average of the smoothness) is initially assigned with p =1, and the region with a smaller smoothness difference (less than the statistical average of the smoothness) is initially assigned with p =2, and then the threshold value ε is fixed to re-calculate the global sparse constraint inversion.

[0024] Preferably, in step S5, the content of the adaptive order dynamic tuning includes:

[0025] In the iterative inversion process, a model update amount discriminant function v is introduced to statistically analyze the difference between the density models before and after two iterations, and a reference value δ (typically, one percent can be selected) is given to define the size of the model update difference; if v is greater than the reference value δ , it is determined that the order p selected in this iteration is still not the optimal value, and needs to be tuned, and if v is less than the reference value δ , the value of p does not need to be adjusted in the next iteration process, and the model update amount discriminant function v is as shown in formula 4:

[0026] Formula 4: ;

[0027] Wherein, m (k) is the density model result of the n th k iteration inversion, m (0) is the initial density model, v (k) represents the model update amount obtained by the n th k iteration inversion calculation, and if v is less than the reference value δ , the tuning is continued, and otherwise, the value of p is not tuned in the subsequent iteration;

[0028] Since a smaller value of p will promote the sparsity of the inversion model solution, and make the model space resolution higher, the tuning method of the value of p is gradually decreased, and the form is shown in formula 5:

[0029] Formula 5: ;

[0030] wherein, p (k) , p (k+1) the first and the first k time iteration k corresponding to the order of the iteration time p .

[0031] The beneficial effects realized by the present application are:

[0032] The present application realizes a sparse constraint inversion method of gravity data by combining model feature driven partitioning and iterative optimization, improves the inversion resolution and accuracy of the geological model, and enhances the applicability of the method technology to complex geological structure interpretation. The density model obtained by the method technology of the present application has improved boundary definition and position accuracy of the blocky geological body, and the range and continuity of the smooth area structure are good. The inversion process reduces the degree of dependence on artificial intervention, can be applied to the actual needs of complex geological structure scenes, and improves the intelligent level of actual gravity data inversion and interpretation. l 2 smooth inversion and globally fixed l 1 sparse inversion results, the boundary definition and position accuracy of the blocky geological body are improved, the range and continuity of the smooth area structure are good, the inversion process reduces the degree of dependence on artificial intervention, can be applied to the actual needs of complex geological structure scenes, and improves the intelligent level of actual gravity data inversion and interpretation. BRIEF DESCRIPTION OF DRAWINGS

[0033] Figure 1 is a flowchart of the model feature driven adaptive sparse constraint inversion method of gravity data.

[0034] Figure 2 is a combination of geological body model space diagram and observed gravity field data.

[0035] Figure 3 is the inversion method of the present application and the conventional fixed p value sparse constraint inversion method comparison effect diagram; in the figure, (a), (b), (c) are respectively the inversion results of the conventional fixed p=1 sparse constraint method, the conventional fixed p=2 sparse constraint method and the inversion method of the present application at the depth of 250m; (d), (e), (f) are respectively the inversion results of the conventional fixed p=1 sparse constraint method, the conventional fixed p=2 sparse constraint method and the inversion method of the present application at y=500m. DETAILED DESCRIPTION

[0036] With reference to the drawings of the embodiments of the present application, the technical solutions in the embodiments of the present application will be described clearly and completely. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments of the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work are within the protection scope of the present application.

[0037] Please refer to Figures 1 to 2 The embodiment of the present application provides a model feature driven gravity data adaptive sparse constraint inversion method, which comprises the following steps:

[0038] S1, grid processing is performed on the collected gravity field data to generate regular grid data and complete three-dimensional discrete grid dissection of underground space;

[0039] The collected gravity field data needs to be gridded, preprocessed into regular grid data with fixed point and line distance arrangement, and then three-dimensional discrete grid dissection of underground space is completed through the horizontal position information of the observation points, which is used for three-dimensional density forward and inversion calculation of the gravity field data.

[0040] S2, a three-dimensional gravity sparse constraint inversion objective function is constructed;

[0041] The three-dimensional gravity sparse constraint inversion objective function is as follows:

[0042] Formula 1: ;

[0043] Wherein, is a gravity field data fitting term, which is used to measure the fitting degree of the observed gravity field data and the gravity field data predicted by the inversion obtained density model, is a model constraint term, which is used to apply prior information constraint conditions to the density model to obtain a reasonable solution in the model space, represents the square operation of the vector l 2 norm, represents the lp power operation of the vector p norm, μ is a regularization parameter, p is a norm order, W d is a data weight matrix, which is usually composed of the inverse of data error, W m is a model weighting function, G is a forward kernel matrix of the gravity field data, d represents the gravity field observation data, m is a model parameter representing the underground medium density distribution respectively;

[0044] To solve the non-differentiability at zero, an approximation of the norm is used lp instead of the strict definition of the norm lp to construct the model constraint term φ m An approximation of (m) is used, which is expressed as follows:

[0045] Equation 2: ;

[0046] where, m i represents each discrete grid cell of the subsurface density model, ε is a small threshold value, ranging from 10 -1 ~10 -6 , to avoid the denominator being equal to zero, enhancing numerical stability.

[0047] S3, model local feature recognition, calculate the smoothness difference of the model;

[0048] First, fix the order p =1, give two sets of threshold values ε of different orders of magnitude difference m 1 and m 2, respectively, complete the gravity data sparse constraint inversion calculation, obtain two independent three-dimensional density solution models ε Since the change rate of the model smoothness is strongly dependent on the inherent characteristics of the corresponding region in the true model, the blocky region that retains sharp edge features shows higher sensitivity to the threshold L ; Therefore, by defining the smoothness of the two inversion model spaces L i , realize the differentiation of different characteristic regions with smoothness difference Δ L , the smoothness measurement operator L i is defined as the sum of the second-order derivatives of the model space, Δ L is the smoothness difference of the inversion model under different threshold values ε , which is expressed as follows:

[0049] Equation 3: .

[0050] where, m i represents each discrete grid cell of the subsurface density model, i represents two sets of threshold values ε of different orders of magnitude difference x independent inversion corresponding to the number, y , z , respectively represent the three orthogonal directions in the Cartesian coordinate system.

[0051] S4. Model feature partitioning constraint: Divide the model region according to the difference in smoothness.

[0052] Based on the difference in smoothness Δ L The size indicates the blocky and smooth features of the model space, and regions with large differences in smoothness (greater than the statistical mean of smoothness) are initially assigned a smoothness value. p =1, initially assigning a smoothness value to regions with small smoothness differences (less than the statistical mean smoothness). p =2, then fix the threshold ε Perform a global sparse constraint inversion calculation again.

[0053] S5, Adaptive order dynamic tuning;

[0054] A discriminant function for model update is introduced during the iterative inversion process. v To statistically analyze the difference between the density models in two iterations, a reference value is given. δ (Typically, one percent is an option) to define the magnitude of the model update discrepancy; if v Greater than the reference value δ Then the order selected for this iteration is determined. p This is still not the optimal value and needs to be optimized. v Less than the reference value δ Then no adjustment is needed in the next iteration. p Value, the discriminant function for model update. v As shown in Formula 4:

[0055] Formula 4: ;

[0056] in, m (k) For the first k The density model results of the next iteration inversion m (0) For the initial density model, v (k) Indicates the first k The model update amount obtained from the next iteration of inversion calculation, if v Less than the reference value δ If the optimization is successful, continue tuning; otherwise, stop tuning in subsequent iterations. p value;

[0057] Due to smaller p The value will promote the sparsity of the solution of the inversion model, resulting in a higher spatial resolution of the model. p The optimization method for the value is to gradually decrease it, as shown in Formula 5:

[0058] Formula 5: .

[0059] in, p (k) , p (k+1) The first k Second and third k The order corresponding to the +1st iteration time p .

[0060] S6. Iterative inversion solution, repeating the dynamic optimization of the order in S5 until the convergence condition is met, and output the final three-dimensional density model.

[0061] This paper implements adaptive sparse constraint inversion calculation of gravity data based on model feature-driven computation. The optimal solution is sought during the sparse constraint inversion iterative process according to the discrimination method in step S5. p The value is adjusted, and the model structure is gradually refined to avoid overly smoothed or over-focused density models, ultimately resulting in a better performance compared to fixed-order models. p The sparse constraint inversion of values ​​yields more stable and reliable underground three-dimensional density structure characteristics.

[0062] Two models simulating a geological body with sharp edges and smooth density variations are presented. The left model has a uniform density, while the right model shows a density that gradually changes from the center outwards. Model illustrations and corresponding gravity field observation data are shown below. Figure 2 Using a fixed order p The sparse constraint inversion method based on the value and the adaptive sparse constraint inversion method based on model feature driving proposed in this invention are used to calculate the distribution model of underground space density parameters, respectively, and the results are as follows: Figure 3 The inversion results, Figure 3 In Figures (a)-(c), the horizontal sections of each inversion result are taken at a depth of 250 m; in Figures (d)-(f), the vertical sections of each inversion result are taken at y=500 m. The results show that the density distribution model obtained by the method of this invention is superior to that obtained by a fixed method. p Value-constrained inversion can be more accurate and has higher resolution. After calculating the root mean square error with the actual density distribution, the accuracy of the global inversion was improved by 12.6%, indicating that the method of this invention has a high degree of automation and stability, and can provide technical support for practical mineral and oil and gas resource exploration applications.

[0063] The above are merely preferred embodiments of the present invention and do not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made using the present invention specification, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. A gravity data adaptive sparse constraint inversion method based on model feature-driven approach, characterized in that, The method includes: S1. The collected gravity field data is processed into a grid to generate regular grid data and complete the three-dimensional discrete grid subdivision of the underground space. S2. Construct the three-dimensional gravity sparse constraint inversion objective function; S3. Local feature identification of the model and calculation of model smoothness differences; S4. Model feature partitioning constraint: Divide the model region according to the difference in smoothness. S5, Adaptive order dynamic tuning; S6. Iterative inversion solution, repeat the dynamic optimization of the order in S5 until the convergence condition is met, and output the final three-dimensional density model. In step S3, the local feature identification of the model and the calculation of the model smoothness difference include: First, through a fixed order p =1, giving two sets of thresholds with different orders of magnitude of difference. ε The sparse constraint inversion calculation of the gravity field data was completed separately to obtain two independent three-dimensional density solution models. m 1 and m 2. By defining the smoothness of the two inversion model spaces L i The difference in smoothness Δ L To achieve differentiation of different feature regions and smoothness L i Defined as the sum of the second derivatives in the model space, Δ L For different thresholds ε The smoothness difference of the inversion model under constraints is expressed in the following form: Formula 3: ； in, m i Each discrete grid cell represents a density model of underground space. i Two sets of thresholds representing differences of different orders of magnitude ε The corresponding number for independent inversion, x , y , z These represent the three orthogonal directions in the Cartesian coordinate system; In step S4, the model feature partitioning constraint, which divides the model regions based on smoothness differences, includes: Based on the difference in smoothness Δ L The size indicates the blockiness and smoothness of the model space, and regions with smoothness differences greater than the statistical mean smoothness are initially assigned a smoothness value. p =1, initially assigning a smoothness value to regions where the smoothness difference is less than the statistical mean smoothness value. p =2, then fix the threshold ε Perform a global sparse constraint inversion calculation again; In step S5, the adaptive order dynamic tuning includes: A discriminant function for model update is introduced during the iterative inversion process. v To statistically analyze the difference between the density models in two iterations, a reference value is given. δ To define the magnitude of the model update difference; if v Greater than the reference value δ Then the order selected for this iteration is determined. p This is still not the optimal value and needs to be optimized. v Less than the reference value δ If so, no adjustment is needed in the next iteration. p Value, the discriminant function for model update. v As shown in Formula 4: Formula 4: ; in, m (k) For the first k The density model results of the next iteration inversion m (0) For the initial density model, v (k) Indicates the first k The model update amount obtained from the next iteration of inversion calculation, if v Less than the reference value δ If the optimization is successful, continue tuning; otherwise, stop tuning in subsequent iterations. p value; p The optimization method for the value is to gradually decrease it, as shown in Formula 5: Formula 5: ; in, p (k) , p (k+1) The first k Second and third k The order corresponding to the +1st iteration time p .

2. The adaptive sparse constraint inversion method for gravity data based on model feature-driven approach according to claim 1, characterized in that, In step S2, the objective function for the three-dimensional gravity sparse constraint inversion takes the following form: Formula 1: ； in, This is a gravity field data fitting term used to measure the degree of fit between the observed gravity field data and the gravity field data predicted by the inverted density model. These are model constraint terms used to impose prior information constraints on the density model in order to obtain a reasonable solution in the model space. Representing vectors l Squaring the 2-norm Representing vectors lp norm p Exponentiation. μ It is a regularization parameter. p For the norm order, W d For the data weight matrix, W m Weighting function for the model, G It is the forward kernel matrix of the gravity field data. d Represents gravitational field observation data, m These are model parameters that characterize the density distribution of underground media; use lp Approximate form of norm lp Norms are rigorously defined and used to construct model constraint terms. φ m The approximate form of (m) is as follows: Formula 2: ; in, m i Each discrete grid cell represents a density model of underground space. ε It is a threshold, ranging from 10. -1 ~10 -6 This is to avoid the denominator being equal to zero, thus enhancing numerical stability.

Images