Gravity data adaptive sparse constraint inversion method based on model feature driving

Abstract

The invention 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, and the method comprises the steps: carrying out the grid processing of collected gravity field data, generating regular grid data, and completing the underground space three-dimensional discrete grid subdivision; constructing a three-dimensional gravity sparse constraint inversion objective function; identifying local features of the model, and calculating the smoothness difference of the model; carrying out model feature partition constraint, and dividing a model region according to smoothness difference; dynamically adjusting and optimizing the self-adaptive order; and carrying out iterative inversion solution, repeatedly executing the order dynamic tuning until a convergence condition is met, and outputting a final three-dimensional density model. According to the method, the sparse constraint inversion method of the gravity data is realized by combining two aspects of model feature driven partitioning and iterative tuning, the inversion resolution and precision of the geologic model are improved, and the applicability of the method technology to explanation of a complex geologic structure is enhanced.

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 gravity data adaptive sparse constraint inversion method based on model feature driving to solve the technical problems mentioned in the background art.

[0005] The gravity data adaptive sparse constraint inversion method based on model feature driving comprises the following steps: 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; S2, a three-dimensional gravity sparse constraint inversion objective function is constructed; S3, model local feature recognition is performed, and the smoothness difference of the model is calculated; S4, model feature partition constraint is performed, and the model area is divided according to the smoothness difference; S5, adaptive order dynamic optimization is performed; S6, iterative inversion solving is performed, the order dynamic optimization of S5 is repeatedly executed until the convergence condition is met, and a final three-dimensional density model is output.

[0006] Preferably, in step S2, the three-dimensional gravity sparse constraint inversion objective function has the following form: Formula 1: ; 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 density model obtained by inversion, 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 , μ is a regularization parameter, p is the 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; In order to solve the zero point non-derivative, the approximate form of lp norm is used instead of lp norm strict definition, and the model constraint term φ m (m) is constructed in the approximate form, which has the following form: Formula 2: ; in, m i Each discrete grid cell representing the underground space density model, ε is a smaller threshold, ranging from 10 -1 ~10 -6 , is to avoid the denominator being equal to zero and enhance numerical stability.

[0007] Preferably, in step S3, the content of identifying local features of the model and calculating the difference in model smoothness includes: First pass the fixed order p =1, giving two sets of thresholds with different orders of magnitude ε , respectively complete the sparse constraint inversion calculation of gravity field data and obtain two independent three-dimensional density solution models m 1 and m 2. Since the rate of change of the model smoothness strongly depends on the inherent characteristics of the corresponding area in the real model, the block area with sharp edge features is more sensitive to the threshold ε shows higher sensitivity; therefore, by defining the smoothness of the two inversion model spaces L i , with the difference in smoothness Δ L , to achieve the distinction between different feature areas, smoothness measurement operator L i Defined as the sum of the second-order derivatives in model space, Δ L For different thresholds ε The smoothness difference of the inversion model under the constraint is as follows: Formula 3: ; in, m i Each discrete grid cell representing the underground space density model, i Two sets of thresholds representing different orders of magnitude differences ε The number corresponding to the independent inversion, x , y , z They represent the three orthogonal directions in the Cartesian coordinate system.

[0008] Preferably, in step S4, the model feature partition constraint, which divides the model area according to the difference in smoothness, includes: According to the smoothness difference Δ L The size of indicates the blockiness and smoothness characteristics of the model space, and the area with large smoothness difference (greater than the statistical smoothness mean) is initially assigned p =1, the area with small smoothness difference (less than the statistical smoothness mean) is initially assigned p= 2, then fix the threshold value ε The global sparse constraint inversion calculation is re-performed.

[0009] Preferably, in step S5, the content of the adaptive order dynamic tuning includes: In the iterative inversion process, a model update amount discriminant function v is introduced to count the difference between the density models before and after iteration, and a reference value δ (typically, one percent can be selected) is given to define the model update difference size; 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 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: Formula 4: ; Wherein, m (k) is the density model result of the nth iteration inversion, k (0) is the initial density model, m (k) , and the model update amount obtained by the nth iteration inversion calculation, if v is less than the reference value k , the tuning is continued, and otherwise, the value is not tuned in the subsequent iteration; v δ Since a smaller value p will promote the sparsity of the inversion model solution and make the model space resolution higher, the tuning method of the value is gradually decreased, which is shown in formula 5: p Formula 5: p ; Wherein, (k) , (k+1) are the orders corresponding to the nth and (n+1)th iterations, respectively. p p k k p

[0010] The beneficial effects realized by the present application are: ​​​​​​The 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 has higher boundary definition and position accuracy of the blocky geological body, better range and continuity of the structure in the smooth area, and reduced dependence on artificial intervention in the inversion process, which can be applied to actual needs of complex geological structure scenes and improves the intelligent level of actual gravity data inversion and interpretation. l 2 Smooth inversion and global fixed l 1 Sparse inversion results, blocky geological body boundary definition and position accuracy are improved, smooth area structure range and continuity are better, and the inversion process reduces the degree of dependence on artificial intervention, which can be applied to complex geological structure scenes and improves the intelligent level of actual gravity data inversion and interpretation. BRIEF DESCRIPTION OF DRAWINGS

[0011] Figure 1 The flowchart of the model feature driven gravity data adaptive sparse constraint inversion method.

[0012] Figure 2 The combination of the geological body model space diagram and the observed gravity field data.

[0013] Figure 3 The inversion method of the 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 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 application at the depth of 250m. DETAILED DESCRIPTION

[0014] The technical solutions in the embodiments of the application will be clearly and completely described below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments of the application. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of the application.

[0015] Please refer to Figures 1 to 2 The embodiment of the application provides a model feature driven gravity data adaptive sparse constraint inversion method, which comprises the following steps: 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; The collected gravity field data needs to be gridded, preprocessed into regular grid data arranged with fixed points and line distance, and then the three-dimensional discrete grid profiling of underground space is completed through the horizontal position information of observation points for three-dimensional density inversion calculation of gravity field data.

[0016] S2, constructing a three-dimensional gravity sparse constraint inversion objective function; The three-dimensional gravity sparse constraint inversion objective function is as follows: Formula 1: ; Wherein, is a gravity field data fitting term, which is used to measure the fitting degree of the observed gravity field data and the predicted gravity field data obtained by inversion of the 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 norm of the vector p , μ is a regularization parameter, p is the 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 gravity field data forward kernel matrix, d represents the gravity field observation data, m is a model parameter representing the density distribution of the underground medium respectively; In order to solve the zero point non-derivative, the approximate form of lp norm is used instead of lp norm strict definition, the model constraint term φ m (m) is constructed in the approximate form, which is expressed as follows: Formula 2: ; Wherein, m i represents each discrete grid element of the underground space density model, ε is a small threshold value, ranging from 10 -1 ~10 -6 , in order to avoid the denominator equal to zero, and enhance the numerical stability.

[0017] S3, model local feature recognition, calculating the smoothness difference of the model; First, by fixing the order p =1, two sets of threshold values of different orders of magnitude difference are givenε , respectively, to obtain two independent 3D density models m 1 and m 2, since the rate of change of the model smoothness strongly depends on the intrinsic characteristics of the corresponding region in the true model, the blocky region with sharp edge features shows higher sensitivity to the threshold ε ; therefore, by defining the smoothness L i , the difference in smoothness L , the smoothness measure operator L i is defined as the sum of the second-order derivatives of the model space, and the difference in smoothness L of the inversion model under different thresholds ε is as follows: Equation 3: .

[0018] wherein m i represents each discrete grid cell of the subsurface density model, i represents two groups of thresholds with different orders of magnitude differences ε , and the corresponding numbers of the independent inversions, x , y , z respectively represent the three orthogonal directions in the Cartesian coordinate system.

[0019] S4, model feature partitioning constraint, the model region is divided according to the difference in smoothness; According to the size of the difference in smoothness L , the blocky and smooth features of the model space are indicated, the region with a larger difference in smoothness (greater than the statistical average of the smoothness) is initially assigned p =1, and the region with a smaller difference in smoothness (less than the statistical average of the smoothness) is initially assigned p =2, and then the threshold ε is fixed to perform sparse constraint inversion calculation in the whole region.

[0020] S5, adaptive order dynamic optimization; In the iterative inversion process, a model update difference function v is introduced to statistically analyze the difference between the density models before and after iteration, and a reference value δ (typically, one percent can be selected) is given to define the model update difference size; 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 optimized, if vLess than the reference value δ , then no adjustment is required in the next iteration p Value, discriminant function of model update amount v As shown in Formula 4: Formula 4: ; in, m (k) For the k The density model results of the iterative inversion are: m (0) is the initial density model, v (k) Indicates the k The model update amount obtained by the iterative inversion calculation is v Less than the reference value δ , then continue to tune, otherwise, no tuning will be done in subsequent iterations p value; Due to the smaller p The value will promote the sparsity of the inversion model solution and make the model spatial resolution higher, so p The value is tuned in a step-by-step manner, as shown in Formula 5: Formula 5: .

[0021] in, p (k) , p (k+1) Respectively k Second and k +1 iteration time corresponding to the order p .

[0022] S6, iterative inversion solution, repeat the order dynamic tuning of S5 until the convergence condition is met, and output the final three-dimensional density model.

[0023] Implement the adaptive sparse constraint inversion calculation of gravity data driven by model characteristics. According to the judgment method of step S5, the optimal p value, and gradually refine the model structure to avoid excessive smoothing or over-focusing of the density model structure. p The sparse constrained inversion of the values ​​can obtain more stable and reliable underground three-dimensional density structure characteristics.

[0024] Simulate a set of two models with sharp edge features and smooth density changes. The density of the geological model on the left is uniform, and the density of the geological model on the right is gradually changing from the center to the surroundings. The model diagram and the corresponding gravity field observation data are as follows: Figure 2 Using a fixed order pThe sparse constraint inversion method based on the value and the adaptive sparse constraint inversion method based on model feature drive proposed in this invention respectively calculate the underground space density parameter distribution model, and obtain the following: Figure 3 The inversion result of Figure 3 (a)-(c) are horizontal slices of the inversion results at a depth of 250m; (d)-(f) are vertical slices of the inversion results at y=500m. The results show that the density distribution model obtained by the inversion method of the present invention is better than that of the fixed p Value-constrained inversion can achieve higher accuracy and resolution. Root mean square error calculations compared to the true density distribution show a 12.6% improvement in global inversion accuracy, demonstrating the method's high degree of automation and stability, providing technical support for practical applications in mineral and oil and gas resource exploration.

[0025] The above are only preferred embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the description of the present invention, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. An adaptive sparse constrained inversion method for gravity data driven by model features, characterized by: The method comprises: S1. Grid the collected gravity field data to generate regular grid data and complete the three-dimensional discrete grid subdivision of the underground space; S2, constructing the three-dimensional gravity sparse constraint inversion objective function; S3, identifying local features of the model and calculating the difference in model smoothness; S4, model feature partition constraint, divides the model area according to the difference in smoothness; S5, dynamic tuning of adaptive order; S6, iterative inversion solution, repeat the order dynamic tuning of S5 until the convergence condition is met, and output the final three-dimensional density model.

2. The method for adaptive sparse constraint inversion of gravity data based on model feature drive according to claim 1 is characterized in that: In step S2, the objective function of the three-dimensional gravity sparse constrained inversion is as follows: Formula 1: ; in, It is the gravity field data fitting item, which is used to measure the fitting degree between the observed gravity field data and the gravity field data predicted by the density model obtained by inversion. is a model constraint term, which is used to impose prior information constraints on the density model to obtain a reasonable solution in the model space. Represents a vector l The square operation of the 2-norm, Represents a vector lp Norm p Power operation, μ is the regularization parameter, p is the norm order, W d is the data weight matrix, W m is the model weighting function, G is the forward kernel matrix of gravity field data, d represents the gravity field observation data, m are model parameters that respectively characterize the density distribution of underground media; use lp The approximate form of the norm is replaced by lp The norm is strictly defined to construct model constraints φ m The approximate form of (m) is expressed as follows: Formula 2: ; in, m i Each discrete grid cell representing the underground space density model, ε is a threshold value ranging from 10 -1 ~10 -6 , is to avoid the denominator being equal to zero and enhance numerical stability.

3. The method for adaptive sparse constraint inversion of gravity data based on model feature drive according to claim 1, characterized in that: In step S3, the content of identifying local features of the model and calculating the difference in model smoothness includes: First pass the fixed order p =1, giving two sets of thresholds with different orders of magnitude ε , respectively complete the sparse constraint inversion calculation of gravity field data and 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 , with the difference in smoothness Δ L , to achieve the distinction of different feature areas, smoothness L i Defined as the sum of the second-order derivatives in model space, Δ L For different thresholds ε The smoothness difference of the inversion model under the constraint is as follows: Formula 3: ; in, m i Each discrete grid cell representing the underground space density model, i Two sets of thresholds representing different orders of magnitude differences ε The number corresponding to the independent inversion, x , y , z They represent the three orthogonal directions in the Cartesian coordinate system.

4. The method for adaptive sparse constrained inversion of gravity data based on model feature drive according to claim 1, characterized in that: In step S4, the model feature partitioning constraint, which divides the model area according to the difference in smoothness, includes: According to the smoothness difference Δ L The size of indicates the blockiness and smoothness characteristics of the model space, and the area with a smoothness difference greater than the statistical smoothness mean is initially assigned p =1, the area with a smoothness difference less than the statistical smoothness mean is initially assigned p =2, then fix the threshold ε Re-perform the sparse constrained inversion calculation for the entire domain.

5. The method for adaptive sparse constraint inversion of gravity data based on model feature drive according to claim 1, characterized in that: In step S5, the content of dynamic tuning of the adaptive order includes: Introducing a discriminant function of the model update amount in the iterative inversion process v , to count the difference between the density models of the two iterations before and after, given a reference value δ To define the size of the model update difference; if v Greater than the reference value δ , then the order selected for this iteration is determined to be p It is still not the optimal value and needs to be tuned. If v Less than the reference value δ , then no adjustment is required in the next iteration p Value, discriminant function of model update amount v As shown in Formula 4: Formula 4: ; in, m (k) For the k The density model results of the iterative inversion are: m (0) is the initial density model, v (k) Indicates the k The model update amount obtained by the iterative inversion calculation is v Less than the reference value δ , then continue to tune, otherwise, no tuning will be done in subsequent iterations p value; p The value is tuned in a step-by-step manner, as shown in Formula 5: Formula 5: ; in, p (k) , p (k+1) Respectively k Second and k +1 iteration time corresponding to the order p .

Images