Adaptive octree three-dimensional magnetic method inversion method fused with abnormal region recognition

Abstract

The invention discloses a self-adaptive octree three-dimensional magnetic method inversion method fused with abnormal region recognition, and belongs to the technical field of three-dimensional magnetic method data inversion processing. Comprising the steps that initial inversion is rapidly completed on a coarse grid, the overall contour of an underground structure is obtained, a fuzzy c-means clustering algorithm is introduced, and a target area with significant magnetic anomaly is automatically extracted from a coarse solution; and constructing a multi-level octree grid in the extracted abnormal region to realize fine subdivision of the complex geological boundary, and further executing high-precision inversion by using a fine model as a new starting point. The adaptive octree grid has the characteristics of high efficiency, small grid number and small data fitting difference. The process of coarse inversion, intelligent identification, local refinement and fine inversion can be circularly executed as required, and computing resources are dynamically allocated. Compared with a traditional global fine grid scheme, the inversion process has the advantages that the inversion precision is maintained and even improved, and meanwhile, the number of grid units is remarkably reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of three-dimensional magnetic data inversion processing, and in particular to an adaptive octree three-dimensional magnetic inversion method integrating abnormal region identification. Background Art

[0002] Magnetic methods, a key tool in geophysical exploration, rely on magnetic field anomalies caused by magnetic differences in the subsurface medium. Magnetic potential field data are observed to infer subsurface structure and the distribution of anomalous bodies. These methods are widely used in fields such as mineral resource exploration and regional geological structure analysis. Magnetic inversion is the core step in the quantitative interpretation of measured magnetic field data, aiming to recover the position, shape, and physical properties of subsurface magnetic bodies. It serves as a crucial bridge between physical observation and geological interpretation. In three-dimensional magnetic inversion, the subsurface is typically discretized into a large number of grid cells, and the physical properties of each grid cell are estimated by minimizing an objective function. However, the subsurface geological bodies are complex in morphology. High-precision delineation of anomalous body boundaries using a regular hexahedral grid requires a detailed meshing of the study area, which rapidly increases the model's degrees of freedom, significantly increases computational complexity, and places higher demands on storage resources. In recent years, numerous studies have focused on improving the computational efficiency of magnetic forward and inversion methods. Changing the model's meshing scheme can fundamentally reduce model size and improve inversion efficiency. The locally refined grid form can use fine grids only in the key areas for meshing, and still use coarse grids for meshing in large areas, which greatly reduces the computing resources occupied by non-critical areas. However, due to the volume difference of the grid, the locally refined grid has unstable convergence during inversion, and the smooth constraint calculation formula of the regular hexahedral grid cannot be used. x 、 y 、 z Applying differential operators between adjacent meshes to construct a smoothing matrix is ​​crucial, so how to incorporate volumetric effects and enforce smoothing constraints is crucial. Furthermore, the question of "where to refine" remains when refining local meshes. Existing methods often focus on extracting key regions in a specific direction (horizontally or vertically), but remain insufficient for full-scale 3D spatial refinement. Therefore, determining the right method to extract key regions and achieve efficient and high-precision 3D magnetic inversion remains a challenging issue. Summary of the Invention

[0003] The purpose of the present invention is to solve the problems of low efficiency and computational waste in existing large-scale three-dimensional magnetic data inversion, and to provide an adaptive octree three-dimensional magnetic inversion method integrating abnormal area identification.

[0004] An adaptive octree 3D magnetic inversion method integrating abnormal region identification includes the following steps: Step 1: Input magnetic observation data, input the initial coarse grid magnetization intensity model, input the number of grid refinements, initial regularization factor, focusing factor, and depth weighting parameter; Step 2: Introduce depth weighting function and smooth focus regularization method to construct the 3D magnetic data inversion objective function ; Step 3: Use the Gauss-Newton method to optimize and solve the 3D magnetic data inversion objective function to obtain the model iteration expression, and perform initial inversion iteration to obtain the initial coarse grid inversion result; Step 4: Combine the fuzzy c-means clustering algorithm and the adaptive extreme value technology to separate the target area and the abnormal area to achieve intelligent extraction of the abnormal area; Step 5: Use the octree grid to refine the abnormal area to obtain a model with local refined grid division; Step 6: For the non-uniformly divided octree grid, smoothing constraints are performed by applying difference operators on adjacent grids respectively; Step 7: Use the locally refined mesh model as the initial model for inversion iteration to obtain high-precision inversion results, and calculate the data fitting error to determine whether to continue the "intelligent identification-local refinement-fine inversion" process.

[0005] In the second step, the inversion objective function of the three-dimensional magnetic data is constructed, and the expression is as follows: ; Among them, the diagonal matrix is the data weighting matrix, is the normalized observation error of the th point; is the forward response; d is the observation data; m is the physical parameter vector; m ref is the reference model; is the regularization factor; is the model weight matrix; Due to the differences in the volume size of the octree grid, the inversion of the refined grid is not very stable. A volume-dependent depth weighting function is used, as shown below: ; in, V j For the j The volume of a grid; r ij For the i Grid to j The distance between observation points; β A constant used to control the strength of the weighting function, usually 0.5 < β <1.5, N is the number of observation data, Mis the number of model mesh divisions; The smooth focus regularization method is used to obtain inversion results with clearer anomaly boundaries while ensuring convergence stability. The expression of the model weighting matrix is ​​as follows: ; ; ; in, is the focusing part of the smooth aggregation regularization; RRR is the smoothing part of the smooth focusing regularization; e is the focusing factor, which is used to control the focusing degree of the inversion model; R x 、R y and R z They are x 、 y 、 z The difference operator applied between adjacent grids in the direction.

[0006] In step 3, the Gauss-Newton method is used to optimize the objective function formula (1), and the following model iteration expression can be obtained: ; Where J is the Jacobian matrix, For the k The amount of change in model parameters for each iteration; In order to improve the efficiency of 3D magnetic inversion, an initial inversion is first performed on a coarse grid to quickly obtain the overall outline of the underground structure, preparing for the subsequent intelligent identification of abnormal areas.

[0007] In step 4, the fuzzy c-means clustering algorithm can automatically calculate the degree of membership of the physical property value of each grid relative to the cluster center based on the cluster center; the adaptive extreme value technology can adaptively determine the initial cluster center and the number of cluster centers; through the fusion of the fuzzy c-means clustering algorithm and the adaptive extreme value technology, the intelligent extraction of abnormal areas can be accurately achieved; The objective function expression of the fuzzy c-means clustering algorithm is as follows: ; ; in, n is the mesh number, p Expressed as the number of clusters, u ij Expressed as i Model parameters m i Relative to the j Cluster centers vj The membership degree, q is the fuzzy factor; Adaptive extreme value technology first selects the value of the area with the most densely distributed physical property values ​​in the model as the background value a, and replaces the physical property values ​​close to a with a to enhance data discrimination; then, it extracts the extreme value points in the model as the center points of the outliers, and uses this to determine the value and number of cluster centers; After determining the value and number of cluster centers, the membership matrix is ​​obtained by optimizing the objective function (7): ; ; in, t ij For the i Model parameters m i Relative to the j Cluster centers v j The Euclidean distance of t ip For the i Physical property values ​​of the grid m i Relative to the p Cluster centers v p The Euclidean distance of Finally, according to the maximum membership principle, the physical property values ​​can be divided into background areas and abnormal areas, realizing the intelligent extraction of abnormal areas in the low-precision inversion model.

[0008] In step 6, since a grid in the octree grid may be adjacent to multiple grids and cannot be calculated using the formula, all grids adjacent to each grid are first found, and then the difference operator is applied to each grid to ensure the smoothness of the inversion result in any direction in space.

[0009] In step 7, the data fitting error is used to quantitatively evaluate the inversion results, and its expression is as follows: ; Among them, d obs is the observed data, d inv is the predicted forward response of the inversion data, err To observe the error, if the inversion result still does not meet the expected data fitting error requirements, the "intelligent identification-local refinement-fine inversion" process will continue to be carried out to dynamically refine the grid and allocate computing resources during the inversion process, thereby gradually improving the inversion accuracy.

[0010] Beneficial effects of the present invention: 1. Based on fuzzy c-means clustering technology and adaptive cluster center determination method, the present invention can accurately distinguish the background area and target area of ​​the model and realize intelligent extraction of abnormal areas.

[0011] 2. In order to solve the inversion instability problem caused by the volume difference of the octree grid, the present invention introduces a depth weighted matrix and a smooth focus regularization method that considers the volume effect, which improves the inversion stability while ensuring the inversion accuracy.

[0012] 3. The adaptive octree grid partitioning process of the present invention can dynamically adjust the range of the extracted abnormal area during the inversion iteration process, gradually refine the inversion grid, and realize the refinement of the three-dimensional grid space. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 is a flow chart of an embodiment of the present invention; Figure 2 A schematic diagram of applying smooth constraints to a regular hexahedral mesh as described in the background of the present invention; Figure 3 A schematic diagram of applying smoothness constraints to an octree grid according to an embodiment of the present invention; Figure 4 This is a three-dimensional schematic diagram of a V-shaped model according to an embodiment of the present invention; Figure 5 This is a top view schematic diagram of the distribution of measuring points of the V-shaped model according to an embodiment of the present invention; Figure 6 : This is a diagram showing the correspondence between the real model of the V-shaped model of this embodiment and the magnetic response of the real model to be fitted, wherein a is a schematic diagram of the real model of the embodiment of the present invention, and b is a diagram showing the magnetic response of the real model to be fitted; Figure 7 A corresponding relationship diagram of a slice of the inversion result of a uniform coarse grid at x = 5 km and a predicted magnetic response of the inversion result of the uniform coarse grid of the V-shaped model of an embodiment of the present invention, wherein a is a slice of the inversion result of the uniform coarse grid at x = 5 km, and b is a predicted magnetic response diagram of the inversion result of the uniform coarse grid; Figure 8 : This is a graph showing the corresponding relationship between the inversion result slice of the uniform fine grid at x = 5 km and the predicted magnetic response of the uniform fine grid inversion result of the V-shaped model according to an embodiment of the present invention, wherein a is a slice of the inversion result of the uniform fine grid at x = 5 km, and b is a predicted magnetic response graph of the inversion result of the uniform fine grid; Figure 9 This is a correspondence diagram between the inversion result slice of the adaptive octree grid at x=5km of the V-shaped model of an embodiment of the present invention and the predicted magnetic response of the adaptive octree grid inversion result, where a is the inversion result slice of the adaptive octree grid at x=5km, and b is the predicted magnetic response diagram of the adaptive octree grid inversion result. DETAILED DESCRIPTION

[0014] See also Figure 1 、 Figures 3 to 9 Shown is an embodiment of the present invention.

[0015] An adaptive octree 3D magnetic inversion method integrating abnormal region identification includes the following steps: Step 1: Input magnetic observation data, input the initial coarse grid magnetization intensity model, input the number of grid refinements, initial regularization factor, focusing factor, and depth weighting parameter; Step 2: Introduce depth weighting function and smooth focus regularization method to construct the 3D magnetic data inversion objective function ; Step 3: Use the Gauss-Newton method to optimize and solve the 3D magnetic data inversion objective function to obtain the model iteration expression, and perform initial inversion iteration to obtain the initial coarse grid inversion result; Step 4: Combine the fuzzy c-means clustering algorithm and the adaptive extreme value technology to separate the target area and the abnormal area to achieve intelligent extraction of the abnormal area; Step 5: Use the octree grid to refine the abnormal area to obtain a model with local refined grid division; Step 6: For the non-uniformly divided octree grid, smoothing constraints are performed by applying difference operators on adjacent grids respectively; Step 7: Use the locally refined mesh model as the initial model for inversion iteration to obtain high-precision inversion results, and calculate the data fitting error to determine whether to continue the "intelligent identification-local refinement-fine inversion" process.

[0016] In the second step, the inversion objective function of the three-dimensional magnetic data is constructed, and the expression is as follows: ; Among them, the diagonal matrix is the data weighting matrix, is the normalized observation error of the th point; is the forward response; d is the observation data; m is the physical parameter vector; m ref is the reference model; α is the regularization factor; is the model weight matrix; Due to the differences in the volume size of the octree grid, the inversion of the refined grid is not very stable. A volume-dependent depth weighting function is used, as shown below: ; in, V j For the jThe volume of a grid; r ij For the i Grid to j The distance between observation points; β A constant used to control the strength of the weighting function, usually 0.5 < β <1.5, N is the number of observation data, M is the number of model mesh divisions; The smooth focus regularization method is used to obtain inversion results with clearer anomaly boundaries while ensuring convergence stability. The expression of the model weighting matrix is ​​as follows: ; ; ; Among them, R e is the focusing part of the smooth aggregation regularization; RRR is the smoothing part of the smooth focusing regularization; e is the focusing factor, which is used to control the focusing degree of the inversion model; R x 、R y and R z They are x 、 y 、 z The difference operator applied between adjacent grids in the direction.

[0017] In step 3, the Gauss-Newton method is used to optimize the objective function formula (1), and the following model iteration expression can be obtained: ; Where J is the Jacobian matrix, For the k The amount of change in model parameters for each iteration; In order to improve the efficiency of 3D magnetic inversion, an initial inversion is first performed on a coarse grid to quickly obtain the overall outline of the underground structure, preparing for the subsequent intelligent identification of abnormal areas.

[0018] In step 4, the fuzzy c-means clustering algorithm can automatically calculate the degree of membership of the physical property value of each grid relative to the cluster center based on the cluster center; the adaptive extreme value technology can adaptively determine the initial cluster center and the number of cluster centers; through the fusion of the fuzzy c-means clustering algorithm and the adaptive extreme value technology, the intelligent extraction of abnormal areas can be accurately achieved; The objective function expression of the fuzzy c-means clustering algorithm is as follows: ; ; in, n is the mesh number, p Expressed as the number of clusters, u ij Expressed as i Model parameters m i Relative to the j Cluster centers v j The membership degree, q is the fuzzy factor; Adaptive extreme value technology first selects the value of the area with the most densely distributed physical property values ​​in the model as the background value a, and replaces the physical property values ​​close to a with a to enhance data discrimination; then, it extracts the extreme value points in the model as the center points of the outliers, and uses this to determine the value and number of cluster centers; After determining the value and number of cluster centers, the membership matrix is ​​obtained by optimizing the objective function (7): ; ; in, t ij For the i Model parameters m i Relative to the j Cluster centers v j The Euclidean distance of t ip For the i Physical property values ​​of the grid m i Relative to the p Cluster centers v p The Euclidean distance of Finally, according to the maximum membership principle, the physical property values ​​can be divided into background areas and abnormal areas, realizing the intelligent extraction of abnormal areas in the low-precision inversion model.

[0019] In step 6, since a grid in the octree grid may be adjacent to multiple grids and cannot be calculated using the formula, all grids adjacent to each grid are first found, and then the difference operator is applied to each grid to ensure the smoothness of the inversion result in any direction in space.

[0020] In step 7, the data fitting error is used to quantitatively evaluate the inversion results, and its expression is as follows: ; Among them, d obs is the observed data, d invis the predicted forward response of the inversion data, err To observe the error, if the inversion result still does not meet the expected data fitting error requirements, the "intelligent identification-local refinement-fine inversion" process will continue to be carried out to dynamically refine the grid and allocate computing resources during the inversion process, thereby gradually improving the inversion accuracy.

[0021] A grid in the octree grid may be adjacent to multiple grids, which cannot be calculated by the formula. Therefore, it is necessary to first find all the grids adjacent to each grid, and then apply the difference operator to each grid to ensure the smoothness of the inversion result in any direction in space, such as Figure 3 shown.

[0022] The adaptive octree three-dimensional magnetic inversion method for integrating abnormal area identification provided by the present invention is verified below.

[0023] Use Figure 4 The V-shaped model shown in the figure is used for model calculation. The background magnetization intensity is 0 A / m, and the study area of ​​10km×10km×6km is divided into 50×50×60 hexahedral units, and the magnetic inclination and declination are set to 70° and 60° respectively. It contains a V-shaped model with a magnetization intensity of 1 A / m. x The directional distribution length is 1.6km (the length of the anomaly is 1.6km, and the x-axis coordinate in the coordinate system is between 4.2km and 5.8km). y The directional distribution length is 3.6km (that is, the length of the anomaly body is 3.6km, and the y-axis coordinate in the coordinate system is between 3.2km and 6.8km), and the top burial depth is about 0.4km. Figure 5 The distribution of observation points is shown from a bird's-eye view, where the 20×20 black dots represent observation points located on the surface.

[0024] Taking the uniform half-space model (0A / m) as the initial model, three-dimensional inversion experiments were carried out using uniform coarse grid, uniform fine grid and adaptive octree grid to verify the effectiveness of adaptive octree inversion. The uniform coarse grid is divided into 20×20×15 hexahedral elements, and the depth weighting parameter is 1. β =1, focusing factor e =0.02, the number of inversion iterations is 30. The uniform fine grid is divided into 50×50×40 hexahedral elements, and the depth weighting parameter β =1.2, focusing factor e =0.03, and the number of inversion iterations is 20. The adaptive octree grid is refined twice, and its initial grid is the same as the uniform coarse grid. In order to ensure that the results obtained by the initial grid inversion can more macroscopically reflect the distribution range of the anomaly and prevent the extracted anomaly area from failing to cover the distribution range of the anomaly, the initial depth weighting parameterβ =0.8, focusing factor e =0.05, the initial iteration number is 10 times; the depth weighting parameter after refinement β =1, focusing factor e =0.02, the number of mesh iterations after each refinement is 5 more than the number of iterations of the previous mesh.

[0025] Figure 6 is a schematic diagram of the real model and the corresponding magnetic response diagram to be fitted. Figures 7 to 9 The inversion model slices for each grid at x = 5 km are shown, with the anomaly range indicated by white lines. The predicted magnetic response maps corresponding to the inversion models for each grid are shown on the right. Here, x and y represent the horizontal directions, z the vertical direction, κ represents the magnetization intensity, and ΔT represents the magnetic anomaly response value.

[0026] Figure 6 a is a schematic diagram of a true model of a V-shaped model according to an embodiment of the present invention, and b is a magnetic response diagram of the true model to be fitted, and the two correspond to each other; Figure 7 a is a slice diagram of the inversion result of the uniform coarse grid at x=5km of the V-shaped model of an embodiment of the present invention, and b is a predicted magnetic response diagram of the inversion result of the uniform coarse grid, and the two correspond to each other; Figure 8 a is a slice diagram of the inversion result of the uniform fine grid at x=5km of the V-shaped model according to an embodiment of the present invention, and b is a predicted magnetic response diagram of the inversion result of the uniform fine grid, and the two correspond to each other; Figure 9 In the figure, a is a slice diagram of the inversion result of the adaptive octree grid at x=5km of the V-shaped model according to an embodiment of the present invention, and b is a predicted magnetic response diagram of the inversion result of the adaptive octree grid, and the two correspond to each other.

[0027] from Figures 6 to 9 It can be seen that the uniform coarse grid cannot well determine the anomaly's extent, while the uniform fine grid inverts the anomaly downwards. The adaptive octree grid more accurately delineates the anomaly's extent. Compared to the uniform fine grid, the adaptive octree grid (computation time: 173.8s; final inversion grid number: 9192; data fit error: 0.594) is much more efficient than the uniform fine grid (computation time: 37291.1s; final inversion grid number: 100000; data fit error: 0.693). It is easy to see that the adaptive octree grid is characterized by high efficiency, a small number of grids, and a small data fit error.

Claims

1. An adaptive octree 3D magnetic inversion method integrating abnormal region identification, characterized in that: The steps include: Step 1: Input magnetic observation data, input the initial coarse grid magnetization intensity model, input the number of grid refinements, initial regularization factor, focusing factor, and depth weighting parameter; Step 2: Introduce depth weighting function and smooth focus regularization method to construct the 3D magnetic data inversion objective function ; Step 3: Use the Gauss-Newton method to optimize and solve the 3D magnetic data inversion objective function to obtain the model iteration expression, and perform initial inversion iteration to obtain the initial coarse grid inversion result; Step 4: Combine the fuzzy c-means clustering algorithm and the adaptive extreme value technology to separate the target area and the abnormal area to achieve intelligent extraction of the abnormal area; Step 5: Use the octree grid to refine the abnormal area to obtain a model with local refined grid division; Step 6: For the non-uniformly divided octree grid, smoothing constraints are performed by applying difference operators on adjacent grids respectively; Step 7: Use the locally refined mesh model as the initial model for inversion iteration to obtain high-precision inversion results. Calculate the data fitting error to determine whether to continue the "intelligent identification - local refinement - fine inversion" process.

2. The adaptive octree 3D magnetic inversion method integrating abnormal region identification according to claim 1 is characterized by: In the second step, the inversion objective function of the three-dimensional magnetic data is constructed, and the expression is as follows: ； Among them, the diagonal matrix is the data weighting matrix, For the i Normalized observation error of each point; is the forward response; d is the observation data; m is the physical parameter vector; m ref is the reference model; α is the regularization factor; is the model weighting matrix; due to the differences in the volume size of the octree grid, the inversion of the refined grid is not very stable, so a volume-related depth weighting function is used, as shown below: ； in, V j For the j The volume of the grid; r ij For the i Grid to j The distance between observation points; β Is a constant used to control the strength of the weighting function, usually 0.5 < β <1.5, N is the number of observation data, M is the number of model mesh divisions; The smooth focus regularization method is used to obtain inversion results with clearer anomaly boundaries while ensuring convergence stability. The expression of the model weighting matrix is ​​as follows: ； ； ； in, is the focusing part of the smooth aggregation regularization; RRR is the smoothing part of the smooth focusing regularization; e is the focusing factor, which is used to control the focusing degree of the inversion model; R x 、R y 、R z They are x 、 y 、 z The difference operator applied between adjacent grids in the direction.

3. The adaptive octree 3D magnetic inversion method integrating abnormal region identification according to claim 1 is characterized by: In step 3, the Gauss-Newton method is used to optimize the objective function formula (1), and the following model iteration expression can be obtained: ； Where J is the Jacobian matrix, For the k The amount of change in model parameters per iteration; In order to improve the efficiency of magnetic three-dimensional inversion, an initial inversion is first performed on a coarse grid to quickly obtain the overall outline of the underground structure, preparing for the subsequent intelligent identification of abnormal areas.

4. The adaptive octree 3D magnetic inversion method integrating abnormal region identification according to claim 1 is characterized by: In step 4, the fuzzy c-means clustering algorithm can automatically calculate the degree of membership of the physical property value of each grid relative to the cluster center based on the cluster center; the adaptive extreme value technology can adaptively determine the initial cluster center and the number of cluster centers; through the fusion of the fuzzy c-means clustering algorithm and the adaptive extreme value technology, the intelligent extraction of abnormal areas can be accurately achieved; The objective function expression of the fuzzy c-means clustering algorithm is as follows: ； ； in, n is the mesh number, p It is expressed as the number of clusters, u ij Expressed as i Model parameters m i Relative to the j Cluster centers v j The membership degree of q is the fuzzy factor; Adaptive extreme value technology first selects the value of the area with the most densely distributed physical property values ​​in the model as the background value a, and replaces the physical property values ​​close to a with a to enhance data discrimination; then, it extracts the extreme value points in the model as the center points of the outliers, and uses this to determine the value and number of cluster centers; After determining the value and number of cluster centers, the membership matrix is ​​obtained by optimizing the objective function (7): ； ； in, t ij For the i Model parameters m i Relative to the j Cluster centers v j The Euclidean distance of t ip For the i Physical property values ​​of the grid m i Relative to the p Cluster centers v p The Euclidean distance of Finally, according to the maximum membership principle, the physical property values ​​can be divided into background areas and abnormal areas, realizing the intelligent extraction of abnormal areas in the low-precision inversion model.

5. The adaptive octree 3D magnetic inversion method integrating abnormal region identification according to claim 1 is characterized by: In step 6, since a grid in the octree grid may be adjacent to multiple grids and cannot be calculated using the formula, all grids adjacent to each grid are first found, and then the difference operator is applied to each grid to ensure the smoothness of the inversion result in any direction in space.

6. The adaptive octree 3D magnetic inversion method integrating abnormal region identification according to claim 1 is characterized by: In step 7, the data fitting error is used to quantitatively evaluate the inversion results, and its expression is as follows: ； Among them, d obs is the observed data, d inv is the predicted forward response of the inversion data, err To observe the error, if the inversion result still does not meet the expected data fitting error requirements, the "intelligent identification-local refinement-fine inversion" process will continue to be carried out to dynamically refine the grid and allocate computing resources during the inversion process, gradually improving the inversion accuracy.

Images