A three-dimensional magnetic vector inversion method based on a spherical coordinate system and a joint modulus constraint of a curl and a gradient tensor

Abstract

The application relates to a spherical coordinate system three-dimensional magnetic vector inversion method based on a combined modulus constraint of a curl and a gradient tensor, which comprises the following steps: acquiring observation magnetic anomaly data of a region to be inverted, and performing spherical prismatic grid dissection on the region to be inverted to construct a sensitivity matrix of each grid unit; constructing a magnetic vector inversion objective function based on the sensitivity matrix; in the construction process, a global discrete gradient operator is determined through a finite difference scale diagonal matrix, and is introduced into a model regularization functional in the magnetic vector inversion objective function, and a gradient tensor combined modulus reweighting matrix is introduced; an independent curl constraint term of the magnetic vector inversion objective function is constructed, a three-dimensional full-component curl matrix in a discrete condition is introduced into the model regularization functional, a curl reweighting matrix based on a curl modulus is adopted, and a comprehensive objective function is formed; an optimization iterative algorithm is adopted to solve the comprehensive objective function, and the size and spatial magnetization direction of the magnetization intensity are output.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a three-dimensional magnetic vector inversion method in spherical coordinates based on the joint modulus constraint of curl and gradient tensor. Background Technology

[0002] Magnetic inversion is an important means of inferring the spatial distribution and physical properties of underground field sources. When processing geophysical anomaly data at the large regional or global scale, in order to consider the influence of the Earth's curvature, it is necessary to construct a three-dimensional underground model using spherical prisms in spherical coordinates. In addition, due to the influence of strong remanence and demagnetization effects, the magnetization direction of the field source usually deviates from the current geomagnetic field. Traditional scalar inversion usually produces serious deviations. Therefore, a three-dimensional magnetic vector inversion method that considers both the inversion magnetization intensity modulus and direction is needed.

[0003] However, conventional magnetic vector inversion has the following intractable drawbacks when using observed magnetic anomaly data: First, because the vector inversion parameters are three times that of the scalar inversion, the ambiguity of the inversion is greatly increased, which further leads to the problem of blurred inversion boundaries; second, existing penalty constraints are mostly applied independently to each projection component direction, which easily causes the boundary or extension direction of the inverted magnetization modulus structure to deviate from the local coordinate axis and cause coordinate axis polarization of the magnetization direction; third, because there is no joint constraint between the magnetization components, the inverted magnetization direction is prone to non-physical local vortices and disordered deflection, deviating from the macroscopically continuous geological characteristics of the magnetization direction within the same geological body.

[0004] Therefore, it is urgent to introduce a novel comprehensive regularization constraint method into the spherical coordinate system magnetic vector inversion framework. This method needs to simultaneously meet multiple control requirements: first, to reduce the offset of the magnetization modulus structure coordinate axis and the polarization of the magnetization direction caused by the independent gradient constraints of the components; second, to suppress local non-physical vortices and strengthen the joint constraints of each magnetization component; and third, to stably restore the magnetization direction while allowing abrupt changes in the physical property boundaries. Summary of the Invention

[0005] To address the problems of existing technologies, the present invention aims to provide a 3D magnetic vector inversion method in spherical coordinates based on joint modulus constraints of curl and gradient tensors. This method aims to solve the technical defects of existing large-scale magnetic vector inversion methods that rely solely on single-observation magnetic anomaly data, including strong ambiguity of solutions, blurred property boundaries, axial deviation of the magnetization modulus structure and polarization of the magnetization direction due to independent component penalties, and the tendency for local vortices and disordered deflections in the magnetization direction. This method achieves precise characterization of property boundaries and stable reconstruction of spatial magnetization directions without relying on any prior geological models or prior remanent magnetic directions.

[0006] To achieve the above objectives, the present invention provides the following solution: A method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors includes: Acquire the observed magnetic anomaly data of the region to be inverted, and perform spherical prism meshing on the region to be inverted to construct the sensitivity matrix of each mesh cell; The magnetic vector inversion objective function is constructed based on the sensitivity matrix. During the construction process, a global discrete gradient operator is determined by a finite difference scale diagonal matrix and introduced into the model regularization functional of the magnetic vector inversion objective function. At the same time, a gradient tensor joint modulus reweighting matrix is ​​introduced. An independent curl constraint term is constructed for the magnetic vector inversion objective function. At the same time, the three-dimensional full-component curl matrix in the discrete case is introduced into the model regularization functional. A curl reweighting matrix based on the curl modulus is used to form a comprehensive objective function that simultaneously constrains the magnetic vector modulus and direction. An optimization iterative algorithm is used to solve the comprehensive objective function, and the magnitude of the magnetization intensity and the spatial magnetization direction of the underground three-dimensional field source are output.

[0007] Optionally, constructing the sensitivity matrix includes: The underground three-dimensional field source body of the region to be inverted is discretized using a spherical prism in a spherical coordinate system, and observation points and integration source points in the three-dimensional spatial volume domain of the field source body are set. A local rectangular coordinate system is established with the observation point as the origin to determine the magnetic anomaly vector components generated by a single spherical prism at the observation point. In determining these magnetic anomaly vector components, a spatial coordinate transformation matrix is ​​constructed based on the spatial positions of the integration source point and the observation point. This spatial coordinate transformation matrix characterizes the coordinate transformation relationship between the radial, latitudinal, and meridional magnetization components in the local spherical coordinate system of the integration source point and the north, east, and vertically downward components in the local rectangular coordinate system of the observation point. When the observed magnetic anomaly data is the total magnetic anomaly data, the unit direction vector of the background geomagnetic field is determined according to the magnetic declination and magnetic inclination of the background geomagnetic field at the observation point, and the directional projection relationship between the magnetic anomaly vector and the total magnetic anomaly in the local rectangular coordinate system of the observation point is established accordingly. The spatial coordinate transformation matrix is ​​combined with the magnetic anomaly forward modeling kernel function, and when the observed magnetic anomaly data is the total magnetic anomaly data, the directional projection relationship is further combined to form a comprehensive forward modeling kernel function corresponding to the radial, latitudinal, and meridional magnetization components in the local spherical coordinate system of the integration source point, respectively. For each spherical prism mesh element, the integral value of the comprehensive forward modeling kernel function of the total magnetic anomaly in the corresponding mesh element is calculated using three-dimensional Gauss-Legend de Gaulle numerical integration, and the integral value is used as the corresponding forward modeling sensitivity matrix element to construct radial, latitudinal, and meridional forward modeling sensitivity sub-matrices, respectively, and then combined to form a comprehensive forward sensitivity matrix to establish the forward mapping relationship between the subsurface three-dimensional magnetization vector model and the observed magnetic anomaly data: ; in, To observe the column vector of magnetic anomaly data, The three block matrices are for synthesizing the forward sensitivity matrix. , and These correspond to the radial, latitudinal, and meridional forward modeling sensitivity submatrices in spherical coordinates, respectively. This is a column vector of a three-dimensional underground magnetization intensity vector model, with its three component vectors. , and These represent the magnetization components of all subsurface mesh elements in the three orthogonal axes in spherical coordinates.

[0008] Optionally, constructing the magnetic vector inversion objective function includes: ; in, Let be the objective function for magnetic vector inversion. For regularization parameters, It is a data fitting term. Let be the model regularization functional in the continuous spherical space Ω.

[0009] Optionally, the finite difference scale diagonal matrix includes: ; in, A constant factor for controlling the continuity of the structure. The first Spatial partitioning step size of each mesh element along the radial, latitudinal, and longitudinal directions. For the first The geocentric distance between the centers of each grid cell For the first The latitude value of the center of each grid subdivision cell. , , These are diagonal matrixes of scale in the radial, latitudinal, and longitudinal directions, respectively. This represents the total number of discrete grid cells.

[0010] Optionally, introducing the gradient tensor joint modulus reweighting matrix includes: ; in, It is a third-order identity matrix. for Global discrete gradient operator for direction, For spatial difference direction index, For the discrete model, the magnetization three-component vector is... For stable parameters, For gradient sparse norm, For spherical coordinate component indices, It is the gradient modulus reweighting matrix.

[0011] Optionally, an independent curl constraint term is constructed for the objective function of the magnetic vector inversion: ; in, For independent curl constraint terms, These are continuous curl components in three orthogonal directions. These are the relative weighting coefficients for the curl constraint term; For the inverted three-dimensional spatial region; A radial weighting function to overcome the depth attenuation effect, It is a three-dimensional magnetization vector field. For generalized norm parameters, For spherical coordinate differential volume element, It is a vector differential operator.

[0012] Optionally, the continuous curl component includes: ; ; ; in, , , These are the radial, latitudinal, and meridional projection components of the magnetic vector, respectively.

[0013] Optionally, the three-dimensional full-component curl matrix includes: ; in, These are the first-order finite difference operator matrices along the three directions; and These are the radial distance and the latitudinal cosine diagonal matrix, respectively.

[0014] Optionally, the curl reweighting matrix based on curl modulus includes: ; in, For generalized norm parameters, It is a third-order identity matrix. For the discrete model, the magnetization three-component vector is... For stable parameters, It is a three-dimensional full-component curl matrix. For spherical coordinate component indices.

[0015] The beneficial effects of this invention are as follows: 1. Effective suppression of local vortices and disordered deflection: This invention introduces curl constraints in the spherical coordinate regularization framework and uses curl modulus reweighting to suppress local vortices and disordered deflection in the inversion magnetization direction.

[0016] 2. Reduces the axial deviation and magnetization polarization of the modulus structure caused by traditional component-independent penalties: This invention introduces a gradient tensor joint modulus reweighting matrix into the gradient constraint term, so that the spatial gradients of the three projected components are no longer reweighted separately, but are jointly reweighted under a unified three-dimensional gradient modulus, thereby reducing the axial bias and magnetization polarization of the magnetization modulus structure.

[0017] 3. Achieved magnetic vector inversion without prior geological models and prior remanent magnetic directions: This invention combines curl constraints with generalized Lp sparse norm constraints, eliminating the need to construct any prior geological models or add prior remanent magnetic directions. Starting directly from observed anomaly data, it can precisely characterize the physical property boundaries of the field source body and recover the spatial magnetization direction while preserving a strong curl response at boundaries such as fractures or contact zones.

[0018] 4. The geometric distortion caused by the plane projection approximation is reduced by taking into account the influence of the Earth's curvature: The present invention is directly based on the spherical coordinate system that conforms to the geometry of the Earth. When processing large-scale or global-scale magnetic anomaly data, the spatial deformation caused by the plane projection approximation is reduced, which ensures the effectiveness of large-scale potential field processing. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 This invention provides a method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors. Figure 2 This is a schematic diagram of a combined model used to verify the effectiveness of the method of the present invention according to an embodiment of the present invention; Figure 3 This is a model magnetic anomaly diagram of superimposed noise according to an embodiment of the present invention; Figure 4 The images show a horizontal slice comparison of the magnetization intensity and spatial vector direction at a depth of 22.5 km using different methods according to embodiments of the present invention; (a) is the inversion result of the conventional method, and (b) is the inversion result of the method of the present invention. Figure 5 The images show a comparison of the magnitude of magnetization intensity at 115.25°E underground longitude and the direction of the spatial vector using different methods in this invention; (a) is the inversion result of the conventional method, and (b) is the inversion result of the method of this invention. Figure 6 The images show a comparison of the equal-area projected polar coordinate kernel density distribution of the magnetization vector direction in the field source grid cell obtained by different methods in this embodiment of the invention; (a) is the inversion result of the conventional method, and (b) is the inversion result of the method of this invention. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0023] like Figure 1 As shown, this embodiment discloses a three-dimensional magnetic vector inversion method in spherical coordinates based on joint modulus constraints of curl and gradient tensors, including: Step S1: Obtain the observed magnetic anomaly data of the area to be inverted, perform spherical prism meshing on the study area (area to be inverted), and construct the sensitivity matrix of the corresponding mesh.

[0024] Step S2: Set the initial model for magnetic vector inversion and construct the basic objective function for magnetic vector inversion based on the sensitivity matrix. The basic objective function includes a data fitting term and a model regularization term. The initial model refers to the initial values ​​of the parameter vector of the model to be inverted when the magnetic vector inversion begins iteration. The model parameter vector consists of the magnetization intensity components in the three orthogonal directions of each grid cell. In this embodiment, the initial value of the model parameter vector is set to zero, so that the magnetization intensity components in the three directions of each grid cell are iteratively updated from zero. The zero initial model is only used to determine the initial iteration point for the inversion calculation.

[0025] Step S3: Construct first-order finite difference operators and corresponding scale diagonal matrices in three directions of the spherical coordinate system, and combine the two to form a global discrete gradient operator; introduce the gradient tensor joint modulus reweighting matrix, and construct independent curl constraint terms, three-dimensional full-component curl matrix and curl modulus reweighting matrix to form a comprehensive objective function.

[0026] Step S4: Solve the comprehensive objective function using an optimization iterative algorithm to output the magnitude of the magnetization intensity and the spatial magnetization direction of the underground three-dimensional field source body in the study area.

[0027] Further, in step S3, the three-dimensional curl constraint term is defined in continuous spherical coordinates as: ; In the formula, These are the relative weighting coefficients for the curl constraint term; For the inverted three-dimensional spatial region; A radial weighting function to overcome the depth attenuation effect, It is a three-dimensional magnetization vector field. For the generalized norm parameter, satisfying ; It is a differential volume element in spherical coordinates.

[0028] Furthermore, the continuous three-dimensional curl operator The orthogonal projection components along the radial, latitudinal, and longitudinal directions are respectively represented as: ; ; ; In the formula, These are the radial, latitudinal, and longitudinal components of the magnetization vector, respectively.

[0029] Furthermore, in step S3, the non-uniformity correction employs a finite-difference scaled diagonal matrix. They are defined as follows: ; in, A constant factor for controlling the continuity of the structure. The first Spatial partitioning step size of each mesh element along the radial, latitudinal, and longitudinal directions. For the first The geocentric distance between the centers of each grid cell For the first The latitude value of the center of each grid subdivision cell. , , These are diagonal matrixes of scale in the radial, latitudinal, and longitudinal directions, respectively. This represents the total number of discrete grid cells.

[0030] Furthermore, the discrete-form three-dimensional full-component curl matrix is ​​constructed. for: ; In the formula, These are the first-order finite difference operator matrices along the three directions; and These are the radial distance and the latitudinal cosine diagonal matrix, respectively.

[0031] Furthermore, in step S3, in order to reduce the axial deviation of the magnetization modulus structure and the polarization of the magnetization direction caused by the independent gradient penalty of each component, the adaptive reweighting matrix includes a gradient tensor joint modulus reweighting matrix. Its definition is: ; In the formula It is a third-order identity matrix. Let ν be the global discrete gradient operator in the direction ν, where ν is the index of the spatial difference direction. For spherical coordinate component indices. For the discrete model, the magnetization three-component vector is... For stable parameters, It is the gradient sparse norm.

[0032] Furthermore, in step S3, the adaptive reweighting matrix also includes a curl modulus reweighting matrix for suppressing non-physical local vortices and disordered deflections. It is based on the discrete curl matrix structure: ; In the formula, For generalized norm parameters.

[0033] In one possible implementation, as a specific computational implementation of the above-described overall method flow, this embodiment further refines the method into the following steps: Step S101: Obtain the observed magnetic anomaly data of the area to be inverted and set the grid partitioning parameters; Step S102: Divide the target area into spherical prism meshes, and calculate the magnetic anomaly vector response of each mesh cell using Gauss-Legendal numerical integration; when the observed data is the total magnetic anomaly, project the magnetic anomaly vector response along the direction of the background geomagnetic field at the observation point, and construct a comprehensive sensitivity matrix. ; Step S103: Set the null model as the initial iterative model and calculate the radial weighting function to overcome depth decay. ; Step S104: To overcome the physical scale inhomogeneity of the spherical mesh as it varies with latitude, longitude, and depth, calculate the diagonal matrix of the mesh physical scale factor. ; Step S105: Based on the scale factor diagonal matrix, construct a global three-dimensional discrete first-order gradient operator through Kronecker product expansion. ; Step S106: Also based on the scale factor diagonal matrix, independently construct a three-dimensional discrete curl constraint operator matrix in spherical coordinates containing orthogonal difference mappings. ; Step S107: Utilize the magnetic vector model of the current iteration Adaptive update of the model's minimum modulus reweighting term matrix And the curl reweighting matrix based on curl modulus, which is the core control mechanism of this invention and the gradient tensor joint modulus reweighting matrix ; Step S108: Combining the aforementioned reweighting matrices, based on the conventional inversion objective function, the discrete curl constraint term is introduced, and the gradient constraint term is reweighted using the gradient tensor joint modulus reweighting matrix to obtain a new comprehensive objective function; Step S109: Optimize the reconstructed objective function using the Gauss-Newton method, calculate the current model update increment, and update the magnetic vector prediction model. ; Step S110: Determine whether the data fit difference of the objective function or the model update amount in the current iteration meets the preset convergence condition. If not, increment the iteration number and return to step S107 to execute the next iteration; if it meets the condition, stop the iteration and proceed to step S111. Step S111: Output the magnetization intensity of each spherical prism mesh element of the final underground three-dimensional field source body, as well as the reconstructed three-dimensional spatial magnetic declination and magnetic inclination distribution.

[0034] The above method is broken down into steps S101 to S111 in the specific execution flow of the computer. To ensure the logical coherence of the mathematical and physical derivation, the above eleven steps are summarized into three core theoretical stages for detailed description. Among them, steps S101 and S102 correspond to the following step one: forward modeling of magnetic vectors in spherical coordinates; steps S103 to S105 correspond to the following step two: the objective function for magnetic vector inversion; and steps S106 to S111 correspond to the following step three: the objective function for inversion with added curl constraints. The detailed implementation is as follows: Step 1: Forward modeling of magnetic vectors in spherical coordinates; In processing large-scale magnetic anomaly data, this invention employs a spherical prism to discretize the underground field source region in a spherical coordinate system to reduce geometric distortion caused by planar projection approximation. (Observation point...) Three-dimensional spatial volume domain of the source body The source point of integration within is .in and The radial distance from the Earth's center. and Latitude and Longitude.

[0035] A local rectangular coordinate system is established with the observation point as the origin, along three orthogonal directions: north (x-axis), east (y-axis), and vertically downward (z-axis). The magnetic anomaly vector component generated by a single spherical prism at the observation point is denoted as... (Component Index) Its specific kernel integral formula is defined as follows: (1); In the formula The permeability of free space, The magnetization vector of the field source, after spatial coordinate transformation, is represented in the local rectangular coordinate system at the observation point. Components along each coordinate axis, where the component index is... , This is a spherical differential volume element. Three-dimensional forward integral kernel function. The specific expression is: (2); In the formula, Kronecker function ( When it is 1, (when it is 0) and Let a and b be the components of the relative position vector between the source point and the observation point in the directions a and b, respectively, where a, b ∈ {x, y, z}; The distance between the source point and the observation point; , , These represent the components of the relative position vector in the x, y, and z directions, respectively. The angle between the center of the sphere and the observation point is denoted as .

[0036] Since the magnetic anomaly calculation in equation (1) requires the use of the magnetization components in the local rectangular coordinate system of the observation point. The underground field source magnetization intensity vector Defined in the local spherical coordinate system of the field source point, therefore a spatial coordinate transformation matrix is ​​introduced. A coordinate transformation is performed to convert the radial, latitudinal, and meridional magnetization components in the local spherical coordinate system of the source point into the north, east, and vertically downward components in the local rectangular coordinate system of the observation point. The mapping relationship is as follows: (3); In the formula, This is the magnetic vector representation of the magnetization vector of the field source projected onto the local Cartesian coordinate system of the observation point; This is a three-component expression of the magnetization intensity of the field source in its own local spherical coordinate system; for The spatial coordinate transformation matrix, whose specific trigonometric function analytical expansion form is defined as: (4); in , Substituting the coordinate transformation relationship shown in equation (3) into equation (1), the spatial coordinate transformation matrix A is used to establish the magnetization component transformation relationship between the local spherical coordinate system of the integration source point and the local rectangular coordinate system of the observation point, converting the radial, latitudinal, and meridional magnetization components at the integration source point into the north, east, and vertically downward components in the local rectangular coordinate system of the observation point. After substituting the coordinate transformation relationship shown in equation (3) into equation (1), the spatial coordinate transformation matrix A and the magnetic anomaly forward modeling integration kernel function... These are combined to form a comprehensive forward kernel function corresponding to the radial, latitudinal, and meridional magnetization components, respectively.

[0037] Furthermore, the underground three-dimensional field source body of the region to be inverted is discretized using a spherical prism in a spherical coordinate system, and observation points and integration source points within the three-dimensional spatial volume domain of the field source body are set. A local rectangular coordinate system is established with the observation point as the origin to determine the magnetic anomaly vector components generated by a single spherical prism at the observation point. In determining the magnetic anomaly vector components, a spatial coordinate transformation matrix is ​​constructed based on the spatial positions of the integration source point and the observation point. The spatial coordinate transformation matrix is ​​used to characterize the coordinate transformation relationship between the radial, latitudinal, and meridional magnetization intensity components in the local spherical coordinate system of the integration source point and the north, east, and vertically downward components in the local rectangular coordinate system of the observation point. When the observed magnetic anomaly data is the total magnetic anomaly data, the unit direction vector of the background geomagnetic field is determined based on the magnetic declination and magnetic inclination of the background geomagnetic field at the observation point, and accordingly... A directional projection relationship is established between the magnetic anomaly vector and the total magnetic anomaly in the local rectangular coordinate system of the observation point. The spatial coordinate transformation matrix is ​​combined with the magnetic anomaly forward modeling kernel function, and when the observed magnetic anomaly data is the total magnetic anomaly data, the directional projection relationship is further combined to form a comprehensive forward modeling kernel function corresponding to the radial, latitudinal, and meridional magnetization intensity components in the local spherical coordinate system of the integration source point, respectively. For each spherical prism grid cell, the integral value of the comprehensive forward modeling kernel function of the total magnetic anomaly in the corresponding grid cell is calculated using three-dimensional Gauss-Legends numerical integration, and the integral value is used as the corresponding forward modeling sensitivity matrix element to construct radial, latitudinal, and meridional forward modeling sensitivity sub-matrices, respectively, and then combined to form a comprehensive forward modeling sensitivity matrix to establish a forward mapping relationship between the underground three-dimensional magnetization intensity vector model and the observed magnetic anomaly data.

[0038] Specifically, from equation (1) we obtain the first... Observation points The three magnetic anomaly components at the location , and The magnetic anomaly vector in the local rectangular coordinate system of the observation point constitutes the magnetic anomaly vector. When the observed data is the total magnetic anomaly, the first The total magnetic anomaly at each observation point is denoted as It is determined by the projection of the magnetic anomaly vector onto the background geomagnetic field direction at the observation point, i.e. ,in, This represents the unit direction vector of the background geomagnetic field at the observation point in the north, east, and vertically downward local rectangular coordinate systems. and These are the magnetic declination and magnetic inclination of the background geomagnetic field at the corresponding observation point, respectively; the magnetic declination is zero in the north and positive in the east, and the magnetic inclination is zero in the horizontal plane and positive downward.

[0039] When performing forward sensitivity calculations, the underground source region is discretized into... Each spherical prism grid cell was used, and the magnetic anomaly response of each grid cell to each observation point was calculated using three-dimensional Gauss-Legends numerical integration. At each numerical integration node, the spatial coordinate transformation matrix A corresponding to the spatial position of the integration source point and the observation point was determined according to equation (4), and it was combined with the magnetic anomaly forward modeling kernel function shown in equation (2). The obtained comprehensive forward modeling kernel function was numerically integrated in each spherical prism grid cell, and the integration result was used as the corresponding forward modeling sensitivity matrix element. Radial, latitudinal, and meridional forward modeling sensitivity sub-matrices were constructed respectively, and combined to form a comprehensive forward modeling sensitivity matrix, resulting in the following forward modeling matrix equation: (5); in To observe the column vector of total magnetic anomaly data, The first one is obtained by projecting the background geomagnetic field direction according to the above-mentioned relationship. The total magnetic anomaly at each observation point This represents the total number of observation points; The three block matrices are for synthesizing the forward sensitivity matrix. , and The dimension is , and respectively correspond to the radial, latitudinal and meridional forward modeling sensitivity submatrices in spherical coordinates; This is a column vector of a three-dimensional underground magnetization intensity vector model, with its three component vectors. , and All dimensions are , and respectively represent the magnetization components of all subsurface mesh elements along three orthogonal axes in spherical coordinates, where This represents the total number of underground field source grid cells in the subdivision.

[0040] Step 2: Objective function for magnetic vector inversion; The magnetic vector inversion in spherical coordinates has strong ambiguity, requiring it to be transformed into a nonlinear optimization problem that minimizes the objective function, which is constructed as follows: (6); In the formula, For regularization parameters, It is the data fitting term, which is defined as follows: (7); In the formula A weighted matrix of data, where For the first The standard deviation of random errors for each observation. To prevent over-smoothing, the Lp norm is introduced in continuous spherical space. The internal model regularization functional is defined as: (8); In the formula, These are the norm parameters that independently control the modulus and gradient sparsity, respectively. Here, represents the relative weight coefficients of each constraint term, specifically the relative weight coefficients of the model magnitude constraint term, radial gradient constraint term, latitudinal gradient constraint term, and meridional gradient constraint term. The radial weighting function to overcome the skin effect is defined as follows: (9); in The height of the observation surface. It is the average radius of the Earth's surface. The radial distance from the center of the grid cell to the Earth's center. This is the radial weighting factor. When discretizing the aforementioned continuous physical quantities, the physical scale of the mesh will undergo non-uniform distortion with changes in latitude and depth. Therefore, this invention introduces a finite-difference scale diagonal matrix. The corrections are defined as follows: (10); In the formula, A constant factor for controlling the continuity of the structure. The total number of discrete grid cells. The first Spatial partitioning step size of each grid subdivision unit along the radial, latitudinal, and longitudinal directions; For the first The geocentric distance between the centers of each grid cell; For the first The latitude values ​​of the centers of each grid subdivision cell. Combining the above scale factors, the first-order finite difference operator in the three directions of the spherical coordinate system. The global discrete gradient operator is obtained by reconstructing the third-order identity matrix using the Kronecker extension. as follows: (11); in, Using spatial difference direction indices, and based on this, the model is solved using iterative reweighted least squares, and the regularized functional equivalent expansion is: (12); In the formula, the comprehensive weight matrix is ​​defined as the product of the spatial geometric weighting matrix and the reweighting matrix: , The spatial geometric weighting matrix , It is a radial weighted vector. For the volume vector of the mesh element, the operator This is the Hadamard product. To mitigate the spatial misalignment of inversion components caused by the traditional component independence penalty, the reweighting matrix of the minimum modulus term in this invention... Using modulus-constrained operator form: (13); in To avoid the stability parameter of division-by-zero singularity, For spherical coordinate component indices. For the gradient reweighting matrix. If we use the Lawson formula for independent decoupling of conventional components as an approximation: (14); This conventional independent gradient penalty, which reweights the spatial gradients of each magnetization component separately, can easily cause the boundaries or extension directions of the inverted magnetization modulus structure to shift towards the local spherical coordinate axes, and cause coordinate axis polarization in the magnetization direction. To address this, this invention improves upon the traditional independent decoupling matrix by introducing a gradient tensor joint modulus reweighting matrix. This matrix uses a unified reweighting control achieved by calculating the combined modulus of the gradient changes in the three spatial directions, specifically defined as: (15); By using the above formula for joint reweighting, the axial deviation of the magnetization modulus structure boundary or extension direction, as well as the polarization of the magnetization direction, can be reduced simultaneously.

[0041] Although the objective function described above has jointly constrained the magnetization modulus and its spatial gradient, it has not explicitly characterized the local rotational changes of the magnetization vector, and the inversion results may still exhibit local vortices and disordered directional deflections. Therefore, it is necessary to further introduce differential structure constraints that can characterize the local rotational changes of the magnetic vector into the regularization framework.

[0042] Step 3: Add curl constraint to the inversion objective function; To suppress local vortices and disordered directional deflections in the inverted magnetization vector, this invention introduces the curl of the magnetization vector as a differential structure regularization term into the model regularization functional of equation (12). For approximately uniform magnetization regions, the spatial variation of the magnetization vector is small, and its curl response is correspondingly weak. Therefore, by penalizing the curl modulus, local vortices and disordered directional deflections in the inversion are suppressed. In spherical coordinates... The continuous curl components in the three orthogonal directions are expressed as follows: (16); (17); (18); Based on the curl components mentioned above, construct the independent curl constraint term for the objective function. for: (19); In the formula To define the relative weighting coefficients of the curl constraint term, a generalized Lp norm parameter is introduced. To characterize the non-uniform abrupt change at the field source boundary. At structural boundaries such as faults or lithological contact zones, the magnetization vector may undergo abrupt changes and produce a strong curl response; therefore, a sparse norm is employed. This allows the curl response to become spatially sparse, while allowing a small number of strong responses at the boundary to be preserved, thereby mitigating excessive smoothing of boundary abrupt changes.

[0043] Based on the continuous differential form of the curl components described above, and combined with the diagonal matrix of the scale factor, a three-dimensional full-component curl matrix in the discrete case is constructed. : (20); After introducing the curl constraint matrix mentioned above, the comprehensive regularization functional of the model is reconstructed as follows: (twenty one); In the formula, the comprehensive weight matrix Curl reweighting matrix The Lawson formula approximation is based on the Lp norm. To suppress non-physical local vortices and disordered deflection, this invention employs a curl reweighting matrix based on the curl modulus. Define it as: (twenty two); Finally, based on the objective function with added curl constraints constructed above, a zero model is used as the initial iterative model, and the Gauss-Newton method is employed for nonlinear iterative solution. During the iteration process, the reweighting matrix of each part is updated using the updated prediction model at each step, and the solution is adaptively corrected until the convergence condition is met, outputting the magnitude and direction of magnetization of the underground three-dimensional field source body.

[0044] Comparison of inversion results based on theoretical models: Theoretical model parameter design: In this embodiment, to verify the ability of the method of the present invention to recover the morphology and magnetization direction of complex remanent magnetization structures, a theoretical model parameter design is constructed as follows: Figure 2 The three-dimensional synthetic theoretical model shown has a horizontal observation area set from 110° to 120° east longitude and from 10° to 22° north latitude, with both the geomagnetic declination and dip angle set to 45°.

[0045] The model comprises two anomalous bodies with distinct magnetization characteristics, one in the south and one in the north, with their actual magnetization set at 1 A / m. The southern anomalous body is a regular, homogeneous block, geographically distributed between 12 and 14 degrees North latitude and at depths ranging from 15 to 35 km below the surface, its magnetization direction consistent with the background geomagnetic field. The northern anomalous body is a stepped structure spanning 17 to 20.5 degrees North latitude, consisting of four radially uplifted plate-like bodies that descend from deep to shallow depths. The northern anomalous body is characterized by significant remanent magnetization, with a magnetization declination of 225° and a magnetic inclination of 45°. This model balances complex geometrical abrupt changes with the non-uniform distribution of magnetization direction.

[0046] The Gauss-Legendal numerical integral was used to calculate the total magnetic anomaly data generated by the model at an observation altitude of 10 km. Gaussian random noise with a standard deviation of 3% of the maximum anomaly amplitude was added to obtain noisy synthetic total magnetic anomaly data, such as... Figure 3 As shown, this data will be used as the input observation data for subsequent inversion.

[0047] Comparative analysis of inversion effects: Without inputting a priori model of the subsurface medium and a priori magnetization direction, and using a zero model as the initial iterative model, inversion was performed using both conventional and the methods proposed in this invention. The differences in performance between the two methods were analyzed using 22.5 km depth slices, 115.25° longitude vertical slices, and equal-area projection polar coordinate kernel density distribution of the magnetization direction.

[0048] A comparison of the magnitude of magnetization and the spatial vector direction at a depth of 22.5 km is shown in the following image. Figure 4 The conventional magnetic vector inversion method in (a) exhibits blurred source boundaries and outward-spreading artifacts, with a maximum reconstructed magnetization of approximately 0.3–0.4 A / m. Furthermore, the reconstructed magnetization directions are relatively disordered, showing local vortices and disordered deflections at the source boundaries, and exhibiting axial deviations towards local spherical coordinate axis polarization. In contrast, the inversion results obtained using the method provided in this invention, such as… Figure 4 As shown in (b), the reconstructed magnetization is close to the actual design value, the geometric boundary is clearer, and the spatial resolution is significantly improved; at the same time, the magnetic vector arrows are more coherent, which significantly reduces non-physical local vortices and axial deviations.

[0049] A comparison of the magnitude of magnetization at underground longitude 115.25°E with the direction of the spatial vector in a perpendicular slice, as shown in the image. Figure 5 The conventional inversion results shown in (a) exhibit obvious longitudinal artifacts in the vertical direction, making it impossible to determine the true underground burial depth; furthermore, the magnetization directions within and at the boundaries are quite chaotic. In contrast, the results obtained by this invention... Figure 5(b) The results shown not only effectively restored the physical property boundary of the model jump, but also the internal magnetization direction was close to the real magnetization direction, thus restoring the model burial depth and magnetization direction well.

[0050] A comparison of the inverted, equal-area projected polar coordinate kernel density distribution with respect to the magnetization vector direction is shown in the following figures: Figure 6 (a) shows the conventional inversion results, where the directional probability distribution is clearly diffuse, with a wide distribution span and an anomalous distribution area in the shape of north-south vertical stripes. The kernel density peaks deviate significantly from the two white asterisks representing the true direction, and the reconstructed direction deviates significantly from the true direction. In contrast, the results obtained by this invention are as follows: Figure 6 As shown in (b), the magnetization directions reconstructed by the inversion are highly concentrated, and their directional probability densities form sharp and independent kernel density peaks near the two true magnetization directions, indicating that the magnetization directions are well recovered. In summary, the above inversion results show that this method can effectively recover the physical properties of subsurface anomalies and stably reconstruct their magnetization directions.

[0051] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors, characterized in that, include: Acquire the observed magnetic anomaly data of the region to be inverted, and perform spherical prism meshing on the region to be inverted to construct the sensitivity matrix of each mesh cell; The magnetic vector inversion objective function is constructed based on the sensitivity matrix. During the construction process, a global discrete gradient operator is determined by a finite difference scale diagonal matrix and introduced into the model regularization functional of the magnetic vector inversion objective function. At the same time, a gradient tensor joint modulus reweighting matrix is ​​introduced. An independent curl constraint term is constructed for the magnetic vector inversion objective function. At the same time, the three-dimensional full-component curl matrix in the discrete case is introduced into the model regularization functional. A curl reweighting matrix based on the curl modulus is used to form a comprehensive objective function that simultaneously constrains the magnetic vector modulus and direction. An optimization iterative algorithm is used to solve the comprehensive objective function, and the magnitude of the magnetization intensity and the spatial magnetization direction of the underground three-dimensional field source are output.

2. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors as described in claim 1, characterized in that, Constructing the sensitivity matrix includes: The underground three-dimensional field source body of the region to be inverted is discretized using a spherical prism in a spherical coordinate system, and observation points and integration source points in the three-dimensional spatial volume domain of the field source body are set. A local rectangular coordinate system is established with the observation point as the origin to determine the magnetic anomaly vector components generated by a single spherical prism at the observation point. In determining these magnetic anomaly vector components, a spatial coordinate transformation matrix is ​​constructed based on the spatial positions of the integration source point and the observation point. This spatial coordinate transformation matrix characterizes the coordinate transformation relationship between the radial, latitudinal, and meridional magnetization components in the local spherical coordinate system of the integration source point and the north, east, and vertically downward components in the local rectangular coordinate system of the observation point. When the observed magnetic anomaly data is the total magnetic anomaly data, the unit direction vector of the background geomagnetic field is determined according to the magnetic declination and magnetic inclination of the background geomagnetic field at the observation point, and the directional projection relationship between the magnetic anomaly vector and the total magnetic anomaly in the local rectangular coordinate system of the observation point is established accordingly. The spatial coordinate transformation matrix is ​​combined with the magnetic anomaly forward modeling kernel function, and when the observed magnetic anomaly data is the total magnetic anomaly data, the directional projection relationship is further combined to form a comprehensive forward modeling kernel function corresponding to the radial, latitudinal, and meridional magnetization components in the local spherical coordinate system of the integration source point, respectively. For each spherical prism mesh element, the integral value of the comprehensive forward modeling kernel function of the total magnetic anomaly in the corresponding mesh element is calculated using three-dimensional Gauss-Legend de Gaulle numerical integration, and the integral value is used as the corresponding forward modeling sensitivity matrix element to construct radial, latitudinal, and meridional forward modeling sensitivity sub-matrices, respectively, and then combined to form a comprehensive forward sensitivity matrix to establish the forward mapping relationship between the subsurface three-dimensional magnetization vector model and the observed magnetic anomaly data: ； in, To observe the column vector of magnetic anomaly data, The three block matrices are for synthesizing the forward sensitivity matrix. , and These correspond to the radial, latitudinal, and meridional forward modeling sensitivity submatrices in spherical coordinates, respectively. This is a column vector of a three-dimensional underground magnetization intensity vector model, with its three component vectors. , and These represent the magnetization components of all subsurface mesh elements in the three orthogonal axes in spherical coordinates.

3. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors as described in claim 1, characterized in that, The objective function for constructing the magnetic vector inversion includes: ； in, Let be the objective function for magnetic vector inversion. For regularization parameters, It is a data fitting term. Let be the model regularization functional in the continuous spherical space Ω.

4. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors according to claim 1, characterized in that, The finite difference scale diagonal matrix includes: ； in, A constant factor for controlling the continuity of the structure. The first Spatial partitioning step size of each mesh element along the radial, latitudinal, and longitudinal directions. For the first The geocentric distance between the centers of each grid cell For the first The latitude value of the center of each grid subdivision cell. , , These are diagonal matrixes of scale in the radial, latitudinal, and longitudinal directions, respectively. This represents the total number of discrete grid cells.

5. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors according to claim 1, characterized in that, The introduction of the gradient tensor joint modulus reweighting matrix includes: ； in, It is a third-order identity matrix. for Global discrete gradient operator for direction, For spatial difference direction index, For the discrete model, the magnetization three-component vector is... For stable parameters, For gradient sparse norm, For spherical coordinate component indices, It is the gradient modulus reweighting matrix.

6. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors according to claim 1, characterized in that, Construct the independent curl constraint term for the magnetic vector inversion objective function: ； in, For independent curl constraint terms, These are continuous curl components in three orthogonal directions. These are the relative weighting coefficients for the curl constraint term; For the inverted three-dimensional spatial region; A radial weighting function to overcome the depth attenuation effect, It is a three-dimensional magnetization vector field. For generalized norm parameters, For spherical coordinate differential volume element, It is a vector differential operator.

7. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors according to claim 6, characterized in that, The continuous curl components include: ； ； ； in, , , These are the radial, latitudinal, and meridional projection components of the magnetic vector, respectively.

8. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors according to claim 1, characterized in that, The three-dimensional full-component curl matrix includes: ； in, These are the first-order finite difference operator matrices along the three directions; and These are the radial distance and the latitudinal cosine diagonal matrix, respectively.

9. The method for three-dimensional magnetic vector inversion in spherical coordinates based on joint modulus constraints of curl and gradient tensors according to claim 1, characterized in that, The curl reweighting matrix based on curl modulus includes: ； in, For generalized norm parameters, It is a third-order identity matrix. For the discrete model, the magnetization three-component vector is... For stable parameters, It is a three-dimensional full-component curl matrix. For spherical coordinate component indices.

Images