A three-dimensional inversion method for airborne magnetic survey based on non-structural node parameterization

Abstract

The application discloses a kind of aviation magnetic method three-dimensional inversion methods based on non-structure node parameterization, this method is in the condition that magnetization direction is fixed and does not consider remanence and demagnetization effect, define magnetization intensity amplitude parameter at non-structure tetrahedron grid node, and construct continuous magnetization model by linear shape function;Equivalent reorganization is node parameterization forward response to unit analytical response, assemble sensitivity matrix by the way of unit contribution to global node column accumulation, construct inversion objective function including node parameterization and model smooth constraint sensitivity matrix, and obtain underground magnetization intensity amplitude distribution by Gauss-Newton and LSQR joint solution.Compared with traditional unit constant parameterization method, the application is helpful to reduce the unknown number scale, improve assembly and solving efficiency, and improve result stability and geological rationality.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, specifically a three-dimensional inversion method based on airborne magnetics with unstructured node parameterization. Background Technology

[0002] Current three-dimensional inversion methods using airborne magnetics generally employ a constant parameterization approach based on grid cells. This involves discretizing the underground space into a large number of regular or unstructured grid cells and assuming that the magnetization amplitude within each cell is a fixed constant. While this modeling approach is relatively straightforward in forward modeling, it introduces two interrelated and difficult-to-reconcile drawbacks in the inversion solution.

[0003] First, the excessively large scale of unknowns in the inversion leads to severe underdeterminacy and strong ambiguity in the inversion problem. To accurately characterize the spatial morphology of complex geological bodies, the number of grid cells is typically in the hundreds of thousands or even millions. Under constant cell parameterization, each cell constitutes an independent parameter to be inverted, making the dimension of the unknowns on the same order of magnitude as the total number of cells. However, the observation data from airborne magnetic surveys are constrained by the spacing between survey lines and flight altitude, resulting in far fewer effective data points than unknowns, forming a severely underdetermined set of inversion equations. This means that even extremely small noise or discretization errors in the observation data can cause the inversion solution to fluctuate drastically within a vast null space. Multiple sets of magnetic distribution models with vastly different characteristics can fit the observation response almost equally, resulting in poor inversion stability and a tendency to generate non-physical spurious anomalies.

[0004] Secondly, the unit constant parameterization mathematically approximates the true magnetic distribution with piecewise constant functions, failing to characterize the continuous variation of underground magnetism and leading to non-physical abrupt changes in the inversion results. The magnetism of actual geological bodies often exhibits a spatially continuous gradual change, such as the decreasing magnetic susceptibility from the center to the edge of an intrusive body or the gradient transition of magnetic minerals within an alteration zone. However, the unit constant model forces the magnetization intensity to remain uniform and constant within the unit, while allowing numerical jumps at the unit boundaries. Essentially, it can only obtain piecewise constant solutions composed of discrete steps. Even with the introduction of smoothing regularization, the results still present a stepped or "block-like" illusion. Artificial boundaries between units are easily misinterpreted as clear geological interfaces, obscuring continuously changing structures such as gradual alteration zones and edge facies zones, deviating from true geological laws. This non-physical abrupt change introduced by the parameterization method itself severely restricts the accuracy and geological rationality of the inversion interpretation. Summary of the Invention

[0005] The purpose of this invention is to provide a three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, the steps of which include:

[0006] S1: Define the magnetization amplitude parameter on the nodes of the unstructured tetrahedral mesh, and interpolate it within the tetrahedral element using a linear shape function to obtain a piecewise linear continuous magnetization model.

[0007] S2: Based on the piecewise linear continuous magnetization model, the analytical response of the linear magnetization of the unit is equivalently recombined to obtain the analytical expression of the forward response consistent with the nodal parameterization.

[0008] S3: Based on the forward response analytical formula, construct the sensitivity matrix by accumulating the local contributions of the units to the global node columns;

[0009] S4: Based on the sensitivity matrix, construct an inversion objective function that includes a data fitting term and a model constraint term based on node parameterization;

[0010] S5: Minimize the norm of the node parameterized model and the norm of the spatial gradient change, solve the constraint terms of the node parameterized model, and thus update the inversion objective function;

[0011] S6: The Gauss-Newton method is used to iteratively linearize the inversion objective function, and an augmented sparse least squares subproblem is constructed by combining the upper and lower bound constraints of the model parameters with variable transformation.

[0012] The LSQR algorithm is used to solve the augmented sparse least squares subproblem to obtain the distribution model of underground magnetization amplitude.

[0013] Furthermore, the piecewise linear continuous magnetization model is shown below:

[0014]

[0015] in, For any spatial location within it The continuous magnetization amplitude at that point, For unit numbering, Numbering of local nodes in the unit. It is a unit The The magnetization amplitude at each local node Representation unit Middle and the first Local scalar linear form functions associated with local nodes; , , , For coefficients;

[0016] Among them, coefficient , , , As shown below:

[0017]

[0018]

[0019]

[0020]

[0021]

[0022] In the formula, Tetrahedral unit Volume; , , and Tetrahedral unit The geometric coordinates of the four vertices; , , , Number the 4 vertices Positive cyclical replacement.

[0023] Furthermore, the steps for equivalent reconstruction of the analytical response to linear magnetization of a single element include:

[0024] observation point Set the origin of the local coordinate system, and any source point in this local coordinate system is represented as... ;

[0025] Construct the magnetization amplitude function, i.e.:

[0026]

[0027] in, , Units The constant term and the linear term of the interpolation function;

[0028] Based on the magnetization amplitude function, a forward-modeling analytical expression consistent with the nodal parameterization is constructed, namely:

[0029]

[0030] in, For nodes The magnetization amplitude parameter is defined at this location. As a constant, Permeability in a vacuum environment , For observation point With polyhedra Any source point The distance between them For the first polyhedron individual aspects Outward unit normal vector Observation point to noodles vertical distance, For the number of faces of a polyhedron, Represent a face The total number of edges on, Representing an edge Outward unit normal vector on, This is the gradient operator.

[0031] Furthermore, the sensitivity matrix In the middle, the first Line number The column element is the first The observation data for the first The first-order partial derivatives of the node parameters, i.e.

[0032]

[0033] in, For nodes When taking the unit amplitude and setting all other nodal parameters to zero, at the observation point The unit magnetic flux density generated at that location;

[0034] When the observed data is total magnetic anomaly data, the elements of the sensitivity matrix are represented as follows:

[0035]

[0036] in For the first The unit vector of the geomagnetic field direction at each observation point.

[0037] Furthermore, the inversion objective function, which includes data fitting terms and node-based parameterization model constraint terms, is shown below:

[0038]

[0039] in, This is the data fit difference term. This represents the norm constraint term of the regularization model based on node parameterization. The regularization factor is used to adjust the weight of the regularization term in the objective function; Represents the magnetization magnitude vector of the grid nodes;

[0040] Furthermore, the data fit difference term is shown below:

[0041]

[0042] in for Diagonal data weight matrix; 1st dimension OK elements Indicates the first Estimate the reciprocal of the noise standard deviation at each observation point; For predicting data vectors; For observation data vectors; The size of the observed data; This is a sensitivity matrix based on node parameterization.

[0043] Furthermore, the steps for solving the model constraint terms based on node parameterization include:

[0044] The model norm constraint term, consisting of the model minimization constraint term and the model spatial gradient minimization constraint term, is expressed as a quadratic stable functional based on the L2 norm, namely:

[0045]

[0046] in, For the minimum constraint term of the model based on the L2 norm, The smoothing constraint term for the L2 norm-based node parameterization model is used to penalize the spatial gradient of the magnetization model;

[0047] Order No. The gradient coefficients of the shape functions corresponding to each local node are: , and The unit volume is Construct the model smoothing constraint matrix as follows:

[0048]

[0049]

[0050]

[0051] in For unit numbering, Representative Unit The Middle The global number of each node;

[0052] Calculate the spatial gradient based on the model smoothness constraint matrix. As shown below:

[0053]

[0054] in, For the node-parameterized model smoothing constraint matrix, For gradient operators;

[0055] Based on spatial gradient The model norm constraint terms are calculated.

[0056] Furthermore, before iteratively linearizing the inverted objective function using the Gauss-Newton method, a depth-weighted matrix is ​​used. With horizontal weighted matrix By weighting the node parameterization model, we obtain:

[0057]

[0058]

[0059] in, For nodes depth, The average height of the observation point. The vertical spatial range of the underground inversion model; The value is the vertical magnetization plane observation data or interpolation data corresponding to the underground grid node;

[0060] Known reference model Introducing regularization constraints, we get:

[0061]

[0062] in, , , and These are the weight coefficients of the model's minimum constraint term and the model's smoothness constraint term, respectively.

[0063] Furthermore, the augmented sparse least squares subproblem is as follows:

[0064]

[0065] in, For the augmented coefficient matrix, To augment the data vector.

[0066] Furthermore, before using the LSQR algorithm to solve the augmented sparse least squares subproblem, the augmented sparse least squares subproblem is linearized.

[0067] The linearized augmented sparse least squares subproblem is shown below:

[0068]

[0069]

[0070] in, , , , These are the upper and lower limits of the model, respectively. This is a diagonal matrix representing the first derivative of the parameter transformation, with diagonal elements of... .

[0071] This invention addresses the problems of traditional airborne magnetics 3D inversion, where the magnetization amplitude parameter is established on grid cells, leading to numerous inversion unknowns, underdeterminacy and multiple solutions in the inverse problem, and the difficulty in representing continuous changes in subsurface magnetic parameters using constant element parameterization. It proposes a parameterization method that establishes the magnetization amplitude parameter on unstructured tetrahedral mesh nodes. Since the number of nodes in an unstructured tetrahedral mesh is much smaller than the number of elements, this parameterization method effectively reduces inversion unknowns and weakens the underdeterminacy and multiple solutions in the inverse problem. Simultaneously, it better adapts to continuous changes in subsurface magnetic parameters, improves computational efficiency, reduces memory consumption, and yields more stable and reasonable inversion results. Attached Figure Description

[0072] Figure 1 This is a flowchart illustrating a three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization, provided as an embodiment of the present invention.

[0073] Figure 2 This is a schematic diagram illustrating the construction of a parameterized continuous model of unstructured tetrahedral mesh nodes provided in an embodiment of the present invention.

[0074] Figure 3 This is a schematic diagram of the inversion algorithm provided in an embodiment of the present invention.

[0075] Figure 4 This is a schematic diagram of a cube model provided in an embodiment of the present invention.

[0076] Figure 5 This is a cross-sectional view of the cube model inversion result provided in an embodiment of the present invention.

[0077] Figure 6 A three-dimensional view of the cube model inversion result provided in an embodiment of the present invention.

[0078] Figure 7 This is a schematic diagram of a Y-shaped vein model provided in an embodiment of the present invention.

[0079] Figure 8 The image shows the inversion results of the Y-type dike model provided in this embodiment of the invention. Detailed Implementation

[0080] The present invention will be further described below with reference to embodiments, but it should not be construed that the scope of the present invention is limited to the following embodiments. Various substitutions and modifications made based on ordinary technical knowledge and common practices in the art without departing from the above-described technical concept of the present invention should be included within the scope of protection of the present invention.

[0081] Example 1:

[0082] A three-dimensional inversion method based on unstructured node parameterization using airborne magnetics, comprising the following steps:

[0083] S1: Define the magnetization amplitude parameter on the nodes of the unstructured tetrahedral mesh, and interpolate it within the tetrahedral element using a linear shape function to obtain a piecewise linear continuous magnetization model.

[0084] S2: Based on the piecewise linear continuous magnetization model, the analytical response of the linear magnetization of the unit is equivalently recombined to obtain the analytical expression of the forward response consistent with the nodal parameterization.

[0085] S3: Based on the forward response analytical formula, construct the sensitivity matrix by accumulating the local contributions of the units to the global node columns;

[0086] S4: Based on the sensitivity matrix, construct an inversion objective function that includes a data fitting term and a model constraint term based on node parameterization;

[0087] S5: Minimize the norm of the node parameterized model and the norm of the spatial gradient change, solve the constraint terms of the node parameterized model, and thus update the inversion objective function;

[0088] S6: The Gauss-Newton method is used to iteratively linearize the inversion objective function, and an augmented sparse least squares subproblem is constructed by combining the upper and lower bound constraints of the model parameters with variable transformation.

[0089] The LSQR algorithm is used to solve the augmented sparse least squares subproblem to obtain the distribution model of underground magnetization amplitude.

[0090] Example 2:

[0091] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as in Example 1, further showing the piecewise linear continuous magnetization model as follows:

[0092]

[0093] in, For any spatial location within it The continuous magnetization amplitude at that point, For unit numbering, Numbering of local nodes in the unit. It is a unit The The magnetization amplitude at each local node Representation unit Middle and the first Local scalar linear form functions associated with local nodes; , , , For coefficients;

[0094] Among them, coefficient , , , As shown below:

[0095]

[0096]

[0097]

[0098]

[0099]

[0100] In the formula, Tetrahedral unit Volume; , , and Tetrahedral unit The geometric coordinates of the four vertices; , , , Number the 4 vertices Positive cyclical replacement.

[0101] Example 3:

[0102] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of Embodiments 1-2, further comprising the step of equivalent reconstruction of the analytical response of linear magnetization of the element, including:

[0103] observation point Set the origin of the local coordinate system, and any source point in this local coordinate system is represented as... ;

[0104] Construct the magnetization amplitude function, i.e.:

[0105]

[0106] in, , Units The constant term and the linear term of the interpolation function;

[0107] Based on the magnetization amplitude function, a forward-modeling analytical expression consistent with the nodal parameterization is constructed, namely:

[0108]

[0109] in, For nodes The magnetization amplitude parameter is defined at this location. As a constant, Permeability in a vacuum environment , For observation point With polyhedra Any source point The distance between them For the first polyhedron individual aspects Outward unit normal vector Observation point to noodles vertical distance, For the number of faces of a polyhedron, Represent a face The total number of edges on, Representing an edge Outward unit normal vector on, This is the gradient operator.

[0110] Example 4:

[0111] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of embodiments 1-3, further including the sensitivity matrix. In the middle, the first Line number The column element is the first The observation data for the first The first-order partial derivatives of the node parameters, i.e.

[0112]

[0113] in, For nodes When taking the unit amplitude and setting all other nodal parameters to zero, at the observation point The unit magnetic flux density generated at that location;

[0114] When the observed data is total magnetic anomaly data, the elements of the sensitivity matrix are represented as follows:

[0115]

[0116] in For the first The unit vector of the geomagnetic field direction at each observation point.

[0117] Example 5:

[0118] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of embodiments 1-4, further comprising the following inversion objective function including data fitting terms and node parameterization-based model constraint terms:

[0119]

[0120] in, This is the data fit difference term. This represents the norm constraint term of the regularization model based on node parameterization. The regularization factor is used to adjust the weight of the regularization term in the objective function; Represents the magnetization magnitude vector of the grid nodes;

[0121] Example 6:

[0122] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of Examples 1-5, further showing the data fitting difference term as follows:

[0123]

[0124] in for Diagonal data weight matrix; 1st dimension OK elements Indicates the first Estimate the reciprocal of the noise standard deviation at each observation point; For predicting data vectors; For observation data vectors; The size of the observed data; This is a sensitivity matrix based on node parameterization.

[0125] Example 7:

[0126] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of embodiments 1-6, further comprising the following steps for solving the model constraint terms based on node parameterization:

[0127] The model norm constraint term, consisting of the model minimization constraint term and the model spatial gradient minimization constraint term, is expressed as a quadratic stable functional based on the L2 norm, namely:

[0128]

[0129] in, For the minimum constraint term of the model based on the L2 norm, The smoothing constraint term for the L2 norm-based node parameterization model is used to penalize the spatial gradient of the magnetization model;

[0130] Order No. The gradient coefficients of the shape functions corresponding to each local node are: , and The unit volume is Construct the model smoothing constraint matrix as follows:

[0131]

[0132]

[0133]

[0134] in For unit numbering, Representative Unit The Middle The global number of each node;

[0135] Calculate the spatial gradient based on the model smoothness constraint matrix. As shown below:

[0136]

[0137] in, For the node-parameterized model smoothing constraint matrix, For gradient operators;

[0138] Based on spatial gradient The model norm constraint terms are calculated.

[0139] Example 8:

[0140] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of embodiments 1-7, further comprising: before iteratively linearizing the inversion objective function using the Gauss-Newton method, utilizing a depth weighting matrix... With horizontal weighted matrix By weighting the node parameterization model, we obtain:

[0141]

[0142]

[0143] in, For nodes depth, The average height of the observation point. The vertical spatial range of the underground inversion model; The value is the vertical magnetization plane observation data or interpolation data corresponding to the underground grid node;

[0144] Known reference model Introducing regularization constraints, we get:

[0145]

[0146] in, , , and These are the weight coefficients of the model's minimum constraint term and the model's smoothness constraint term, respectively.

[0147] Example 9:

[0148] A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, with the same technical content as any one of embodiments 1-8, further comprising the augmented sparse least squares subproblem as shown below.

[0149]

[0150] in, For the augmented coefficient matrix, To augment the data vector.

[0151] Example 10:

[0152] A three-dimensional inversion method based on unstructured node parameterization of airborne magnetics, with the same technical content as any one of embodiments 1-9, further wherein the augmented sparse least squares problem is linearized before solving the augmented sparse least squares problem using the LSQR algorithm;

[0153] The linearized augmented sparse least squares subproblem is shown below:

[0154]

[0155]

[0156] in, , , , These are the upper and lower limits of the model, respectively. This is a diagonal matrix representing the first derivative of the parameter transformation, with diagonal elements of... .

[0157] Example 11:

[0158] A three-dimensional inversion method based on unstructured node parameterization using airborne magnetics, comprising the following steps:

[0159] S1: By defining magnetization amplitude parameters on unstructured tetrahedral mesh nodes and interpolating within the tetrahedral elements using linear shape functions, a piecewise linear continuous magnetization model is obtained.

[0160] S2: Based on the piecewise linear continuous magnetization model, by equivalently recombining the analytical response of the linear magnetization of the unit, the analytical expression of the forward response consistent with the nodal parameterization is obtained;

[0161] S3: Construct a sensitivity matrix based on the forward response analytical expression, and complete the sensitivity matrix assembly by accumulating the local contributions of the units to the global node columns;

[0162] S4: Based on the sensitivity matrix, construct an inversion objective function that includes a data fitting term and a model constraint term based on node parameterization.

[0163] S5: Minimize the norm of the nodal parameterized model and the norm of the spatial gradient change, solve the model constraint terms of the nodal magnetization amplitude parameter, and construct the model smoothing constraint matrix based on nodal parameterization.

[0164] S6: The Gauss-Newton method is used to iteratively linearize the inversion objective function, construct an augmented sparse least squares subproblem that combines upper and lower bound constraints of model parameters with variable transformation, and then the LSQR algorithm is used to solve it to obtain the underground magnetization intensity amplitude distribution model.

[0165] For any tetrahedral element In the magnetization direction Under the condition of constant consistency, by defining magnetization amplitude parameters on the nodes of the unstructured tetrahedral mesh and interpolating within the tetrahedral elements using linear shape functions, a piecewise linear continuous magnetization model is obtained:

[0166]

[0167] in, For any spatial location within it The continuous magnetization amplitude at that point, For unit numbering, Numbering of local nodes in the unit. It is a unit The The magnetization amplitude at each local node Representation unit Middle and the first Local scalar linear functions associated with local nodes.

[0168] Based on the aforementioned piecewise linear continuous magnetization model, by equivalently recombining the analytical responses of the element linear magnetization, a forward-modeling analytical expression consistent with the nodal parameterization is obtained. The steps to rewrite the magnetization amplitude function within the range as the sum of constant and linear terms include: [details of the steps are missing from the provided text]. If the origin is set as the local coordinate system origin, then any source point in that local coordinate system is represented as... The magnetization amplitude function is then expressed as:

[0169]

[0170] in, , Units The constant term and the linear term of the interpolation function.

[0171] node Interpolate the observation points individually within their topological support domain. The total magnetic flux density response can be written as:

[0172]

[0173] in, For nodes The magnetization amplitude parameter is defined at this location. As a constant, Permeability in a vacuum environment , For observation point With polyhedra Any source point The distance between them For the first polyhedron individual aspects Outward unit normal vector Observation point to noodles vertical distance, For the number of faces of a polyhedron, Represent a face The total number of edges on, Representing an edge Outward unit normal vector on, This is the gradient operator.

[0174] The steps for assembling the sensitivity matrix based on the analytical expression of the forward response, and for completing the sensitivity matrix assembly by accumulating the local contributions of each element to the global node columns, include:

[0175] Let the sensitivity matrix be set. The Line number Column element is defined as the first The observation data for the first The first-order partial derivatives of the node parameters, i.e.

[0176]

[0177] in, For nodes When taking the unit amplitude and setting all other nodal parameters to zero, at the observation point The unit magnetic flux density generated at that location. When the observed data is total magnetic anomaly data, the sensitivity matrix elements are represented as...

[0178]

[0179] in For the first The unit vector of the geomagnetic field direction at each observation point. In the actual assembly process, for each observation point, all tetrahedral elements are traversed sequentially. First, the constant part kernel function and the kernel function of the first term in the three coordinate directions of the element at the observation point are calculated. Then, the local contribution is accumulated into the matrix column of the corresponding global node number by combining the shape function coefficients of the four local nodes in the element.

[0180] Based on the sensitivity matrix, an inversion objective function is constructed that includes a data fitting term and a model constraint term based on node parameterization:

[0181]

[0182] in This is the data fit difference term. This represents the norm constraint term of the regularization model based on node parameterization. The regularization factor is used to adjust the weight of the regularization term in the objective function. This represents the magnitude vector of magnetization intensity at the grid nodes.

[0183] The expression for the data fit difference term is as follows:

[0184]

[0185] in for The dimensional diagonal data weight matrix, the first OK elements Indicates the first Estimate the reciprocal of the noise standard deviation at each observation point, where the noise is a certain proportion of Gaussian noise. To predict data vectors, For the observed data vector, For the size of the observed data, This is the sensitivity matrix based on node parameterization.

[0186] The steps to minimize the norm and spatial gradient variation norm of the nodal parameterized model, solve for the model constraint terms of the nodal magnetization magnitude parameters, and construct the model smoothing constraint matrix based on nodal parameterization include:

[0187] The model norm constraint term consists of the model minimization constraint term and the model spatial gradient minimization constraint term, and is expressed as a quadratic stable functional based on the L2 norm:

[0188]

[0189] in, For the minimum constraint term of the model based on the L2 norm, The L2 norm-based node parameterization model uses a smoothing constraint term to penalize the spatial gradient of the magnetization model:

[0190]

[0191] in, For the node-parameterized model smoothing constraint matrix, For the gradient operator. For the... A tetrahedral unit, which in The text occupies three rows, corresponding to the model gradients. , and Three directional components. If the first... The gradient coefficients of the shape functions corresponding to each local node are: , and The unit volume is Then we have:

[0192]

[0193]

[0194]

[0195] in For unit numbering, Representative Unit The Middle The global number of each node.

[0196] Introducing a depth-weighted matrix With horizontal weighted matrix The node parameterization model was weighted, and the construction of the depth weighting function was adjusted to accommodate the change in parameterization method:

[0197]

[0198]

[0199] in, For nodes depth, The average height of the observation point. The vertical spatial range of the underground inversion model. This can avoid numerical instability issues such as the weight function having a denominator of 0 or a negative value at the bottom of the node parameterized model. The values ​​are the vertical magnetization plane observation data or interpolated data corresponding to the underground grid nodes. For magnetic anomaly data under tilted magnetization, polarization processing (RTP) should be performed first before using it to construct the horizontal weighting matrix. Combined weights They are used together for weighting the node parameterization model.

[0200] If the reference model is known This can be introduced into the regularization constraint to mitigate the non-uniqueness of the inversion. The regularization constraint is then expressed as:

[0201]

[0202] in, , , and These are the weight coefficients of the model's minimum constraint term and the model's smoothness constraint term, respectively.

[0203] The steps for iteratively linearizing the inversion objective function using the Gauss-Newton method, constructing an augmented sparse least squares subproblem that combines upper and lower bound constraints on model parameters with variable transformation, and then solving it using the LSQR algorithm to obtain the underground magnetization amplitude distribution model include:

[0204] The objective function is:

[0205]

[0206] Reconstructed into linear least squares form:

[0207]

[0208] in, For the augmented coefficient matrix, To augment the data vector.

[0209] Based on the augmented linear objective function, a logarithmic transformation is introduced to apply upper and lower limits of physical properties.

[0210]

[0211] in, , These are the upper and lower limits of the model, respectively. The objective function is then transformed into a function about unconstrained variables. Nonlinear least squares problem:

[0212]

[0213] The Gauss-Newton method is used to perform local linearization iteration on the nonlinear objective function, and at the th... In the next iteration, In the current model If we perform a first-order Taylor expansion at and ignore higher-order terms, the objective function can be approximated as the following linear function:

[0214]

[0215] in, , , This is a diagonal matrix representing the first derivative of the parametric transformation, whose diagonal elements are calculated using the chain rule. .

[0216] In each Gauss-Newton iteration, the regularization factor is selected by a line search based on the Occam inversion idea. Within a preset search range, the regularization factor for the current iteration is adaptively determined based on the data fitting error, the error increase, or the error decrease rate.

[0217] In each Gauss-Newton iteration, multiple regularization factors are selected based on the Occam inversion idea. The LSQR algorithm is then used sequentially to solve the linearized objective function under different regularization factors. The optimal regularization factor and its corresponding solution are selected as the solution for this Gauss-Newton iteration based on preset conditions. The detailed implementation rules are as follows: In the... In the next iteration, the regularization factor From the initial maximum value according to Decrease and calculate each The corresponding model vector and nRMS. When the preset condition is met, nRMS = When the value drops below 1, nRMS increases twice consecutively, or the rate of decrease of nRMS is below a threshold (e.g., 0.01), the regularization factor corresponding to the minimum nRMS is selected as the regularization factor for this iteration, and the corresponding model vector is retained.

[0218] Example 12:

[0219] Verification of a three-dimensional inversion method based on unstructured node parameterization using airborne magnetics, including:

[0220] This invention provides a method for three-dimensional magnetotelluric inversion of continuously varying conductivity in blocks, as described in the following embodiments: Figure 1 , Figure 2 and Figure 3 As shown, it includes the following steps:

[0221] S1: By defining magnetization amplitude parameters on unstructured tetrahedral mesh nodes and interpolating within the tetrahedral elements using linear shape functions, a piecewise linear continuous magnetization model is obtained.

[0222] In this embodiment of the invention, in the magnetization direction Under the condition of constant consistency, by defining magnetization amplitude parameters on the nodes of the unstructured tetrahedral mesh and interpolating within the tetrahedral elements using linear shape functions, a piecewise linear continuous magnetization model is obtained:

[0223]

[0224] in, For any spatial location within it The continuous magnetization amplitude at that point, For unit numbering, Numbering of local nodes in the unit. It is a unit The The magnetization amplitude at each local node Representation unit Middle and the first Local scalar linear functions associated with local nodes.

[0225] S2: Based on the piecewise linear continuous magnetization model, the analytical expression of the forward response consistent with the nodal parameterization is obtained by equivalently recombining the analytical response of the unit linear magnetization.

[0226] In this embodiment of the invention, based on the piecewise linear continuous magnetization model, by equivalently recombining the analytical responses of the unit linear magnetization, a forward modeling analytical expression consistent with the nodal parameterization is obtained, and the observation points are... If the origin is set as the local coordinate system origin, then any source point in that local coordinate system is represented as... The magnetization amplitude function is then expressed as:

[0227]

[0228] in, , Units The constant term and the linear term of the interpolation function.

[0229] node Interpolate the observation points individually within their topological support domain. The total magnetic flux density response can be written as:

[0230]

[0231] in, For nodes The magnetization amplitude parameter is defined at this location. As a constant, Permeability in a vacuum environment , For observation point With polyhedra Any source point The distance between them For the first polyhedron individual aspects Outward unit normal vector Observation point to noodles vertical distance, For the number of faces of a polyhedron, Represent a face The total number of edges on, Representing an edge Outward unit normal vector on, This is the gradient operator.

[0232] S3: Construct a sensitivity matrix based on the forward response analytical formula, and complete the sensitivity matrix assembly by accumulating the local contributions of the units to the global node columns.

[0233] In this embodiment of the invention, a sensitivity matrix is ​​set. The Line number Column element is defined as the first The observation data for the first The first-order partial derivatives of the node parameters, i.e.:

[0234]

[0235] in, For nodes When taking the unit amplitude and setting all other nodal parameters to zero, at the observation point The unit magnetic flux density generated at that location. When the observed data is total magnetic anomaly data, the sensitivity matrix elements are represented as...

[0236]

[0237] in For the first The unit vector of the geomagnetic field direction at each observation point. In the actual assembly process, for each observation point, all tetrahedral elements are traversed sequentially. First, the constant part kernel function and the kernel function of the first term in the three coordinate directions of the element at the observation point are calculated. Then, the local contribution is accumulated into the matrix column of the corresponding global node number by combining the shape function coefficients of the four local nodes in the element.

[0238] S4: Based on the sensitivity matrix, construct an inversion objective function that includes a data fitting term and a model constraint term based on node parameterization.

[0239] In this embodiment of the invention, the objective function is expressed as:

[0240]

[0241] in This is the data fit difference term. This represents the norm constraint term of the regularization model based on node parameterization. The regularization factor is used to adjust the weight of the regularization term in the objective function. This represents the magnitude vector of magnetization intensity at the grid nodes.

[0242] The expression for the data fit difference term is as follows:

[0243]

[0244] in for The dimensional diagonal data weight matrix, the first OK elements Indicates the first Estimate the reciprocal of the noise standard deviation at each observation point, where the noise is Gaussian noise. To predict data vectors, For the observed data vector, For the size of the observed data, This is the sensitivity matrix based on node parameterization.

[0245] S5: Minimize the norm of the nodal parameterized model and the norm of the spatial gradient change, solve the model constraint terms of the nodal magnetization amplitude parameter, and construct the model smoothing constraint matrix based on nodal parameterization.

[0246] In this embodiment of the invention, the model norm constraint term consists of a model minimization constraint term and a model spatial gradient minimization constraint term, and is represented as a quadratic stable functional based on the L2 norm:

[0247]

[0248] in, For the minimum constraint term of the model based on the L2 norm, The L2 norm-based node parameterization model uses a smoothing constraint term to penalize the spatial gradient of the magnetization model:

[0249]

[0250] in, For the node-parameterized model smoothing constraint matrix, For the gradient operator. For the... A tetrahedral unit, which in The text occupies three rows, corresponding to the model gradients. , and Three directional components. If the first... The gradient coefficients of the shape functions corresponding to each local node are: , and The unit volume is Then we have:

[0251]

[0252]

[0253]

[0254] in For unit numbering, Representative Unit The Middle The global number of each node.

[0255] Introducing a depth-weighted matrix With horizontal weighted matrix The node parameterization model was weighted, and the construction of the depth weighting function was adjusted to accommodate the change in parameterization method:

[0256]

[0257]

[0258] in, For nodes depth, The average height of the observation point. The vertical spatial range of the underground inversion model. This can avoid numerical instability issues such as the weight function having a denominator of 0 or a negative value at the bottom of the node parameterized model. The values ​​are the vertical magnetization plane observation data or interpolated data corresponding to the underground grid nodes. For magnetic anomaly data under tilted magnetization, polarization processing (RTP) should be performed first before using it to construct the horizontal weighting matrix. Combined weights They are used together for weighting the node parameterization model.

[0259] If the reference model is known This can be introduced into the regularization constraint to mitigate the non-uniqueness of the inversion. The regularization constraint is then expressed as:

[0260]

[0261] in, , , and These are the weight coefficients of the model's minimum constraint term and the model's smoothness constraint term, respectively.

[0262] S6: The Gauss-Newton method is used to iteratively linearize the inversion objective function, construct an augmented sparse least squares subproblem that combines upper and lower bound constraints of model parameters with variable transformation, and then the LSQR algorithm is used to solve it to obtain the underground magnetization intensity amplitude distribution model.

[0263] In this embodiment of the invention, the objective function is reconstructed into a linear least squares form:

[0264]

[0265] in, For the augmented coefficient matrix, To augment the data vector, a logarithmic transformation is introduced based on the augmented linear objective function to apply upper and lower limits of physical properties.

[0266]

[0267] The objective function is then transformed into one concerning unconstrained variables. Nonlinear least squares problem:

[0268]

[0269] The Gauss-Newton method is used to perform local linearization iteration on the nonlinear objective function, and at the th... In the next iteration, In the current model If we perform a first-order Taylor expansion at and ignore higher-order terms, the objective function can be approximated as the following linear function:

[0270]

[0271] in, , , This is a diagonal matrix representing the first derivative of the parametric transformation, whose diagonal elements are calculated using the chain rule. .

[0272] In each Gauss-Newton iteration, the regularization factor is selected by a line search method based on the Occam inversion idea. Within the preset search range, the regularization factor for the current iteration is adaptively determined according to the data fitting error increase or the error decrease rate.

[0273] In each Gauss-Newton iteration, multiple regularization factors are selected based on the Occam inversion idea, and the LSQR algorithm is used sequentially to solve the linearized objective function under different regularization factors. The optimal regularization factor and the corresponding solution are selected as the solution for this Gauss-Newton iteration based on preset conditions.

[0274] This application uses a cube model as an example in its embodiments, such as... Figure 4 As shown, the study area of ​​this model is a rectangular region of 1000m × 1000m horizontally and 500m vertically. The magnetic anomaly model is a regular cube with sides of 200m, buried 250m underground at its center in the study area (as shown in the figure). Its induced magnetization is 2A / m, and its magnetization tilt (I) and deflection (D) are 60° and 30°, respectively. The geomagnetic tilt and deflection are set to 90° and 0°, respectively. The magnitude of the geomagnetic field is taken as 50000nT. Remanence and demagnetization effects are not considered in the numerical simulation experiments of this paper. The observation plane covers the entire study area, with an observation height of 1m above the ground surface and a measuring point spacing of 50m, for a total of 21 × 21 = 441 measuring points. The theoretical magnetic anomaly field is obtained by forward modeling the above-mentioned cubic magnetic body on the observation plane, and then the absolute value of the theoretical observation with a standard deviation of 2% plus 0.5nT of Gaussian noise is added to construct noisy observation data. The study area was divided into 21,251 unstructured tetrahedral elements using TetGen meshing software. The depth-weighted index was then used. and horizontal weighted index The magnetization amplitudes of the model nodes were set to 1.5 and 0.5 respectively, and the upper and lower limits of the magnetization amplitude were set to −1 A / m and 10 A / m respectively. The stop tolerance of the LSQR solver was set to a stricter value, atol=10. −10 btol=10 −15 Finally, under the same mesh partitioning, data conditions, and inversion parameter settings, the method of this invention is compared with the traditional element constant parameterization inversion method. The inversion results are as follows: Figure 5 and Figure 6As shown, both parameterization methods recovered the approximate spatial extent of the anomaly. Notably, under the same tilt magnetization conditions, the results recovered using the proposed algorithm better match the true model in terms of the anomaly's center location and overall geometry, and the inversion results of the proposed algorithm are smoother.

[0275] This embodiment uses a Y-shaped dike model as an example, such as Figure 7 As shown, the model consists of two inclined plate-like bodies, arranged within a study area of ​​x=[0,1000] and y=[0,1000], with a vertical study depth of 500m. The left plate-like body extends shorter along the dip direction, with a top burial depth of 50m and a bottom burial depth of 250m. The right plate-like body extends longer along the dip direction, with a top burial depth of 50m and a bottom burial depth of 400m. The same observation system as in Example 1 was used, with the standard deviation added to the absolute value of the theoretical observation. Gaussian noise. The inverted region is divided into 53,138 unstructured tetrahedral elements and 10,755 nodes; the depth-weighted index is used. and horizontal weighted index Set the values ​​to 1.5 and 0.5 respectively, and uniformly set the LSQR stop tolerance to atol=10. −10 btol=10 −15 The parameters of the two methods in the Y-type dike model inversion are shown in Table 1. The node parameterization method inverted 10755 unknowns, with a total inversion time of 71.3064 seconds and a maximum memory usage of 0.3855 GB. The element mean parameterization method inverted 53138 unknowns, with a total inversion time of 737.485728 seconds and a maximum memory usage of 2.08928 GB.

[0276] Table 1. Parameter information for the two methods in the Y-type dike model inversion.

[0277] Number of unknowns Inversion time (s) Maximum memory (GB) Node parameterization method 10755 71.3064 0.3855 Cell mean parameterization method 53138 737.485728 2.08928

[0278] Inversion results as follows Figure 8 As shown, both methods can generally recover the approximate location and spatial extent of anomalies, but the method of this invention achieves clearer anomaly boundaries and a more compact 3D anomaly morphology, with a significant reduction in non-physical spikes. The number of unknowns required by the element mean parameterization method is approximately five times that of the node parameterization method, while the inversion time consumed by the node parameterization method is only about 1 / 10 of that of the element mean parameterization method, resulting in an inversion speed improvement of approximately one order of magnitude. Furthermore, compared to conventional algorithms, the maximum memory required for the algorithm proposed in this paper is significantly reduced.

[0279] The beneficial effects of this invention are as follows: Addressing the problems of traditional airborne magnetics 3D inversion where the magnetization amplitude parameter is established on grid cells, resulting in numerous inversion unknowns, strong underdeterminism and multiple solutions in the inverse problem, and the difficulty in characterizing continuous changes in subsurface magnetic parameters using constant element parameterization, this invention proposes a parameterization method that establishes the magnetization amplitude parameter on unstructured tetrahedral grid nodes. Since the number of nodes in an unstructured tetrahedral grid is much smaller than the number of elements, this parameterization method can effectively reduce inversion unknowns and weaken the underdeterminism and multiple solutions in the inverse problem. Simultaneously, it can better adapt to continuous changes in subsurface magnetic parameters, improve computational efficiency, reduce memory consumption, and obtain more stable and reasonable inversion results.

[0280] The above descriptions are merely embodiments of the present invention. Commonly known technical solutions or characteristics are not described in detail here. It should be noted that those skilled in the art can make various modifications and improvements without departing from the present invention, and these should also be considered within the scope of protection of the present invention. These modifications and improvements will not affect the effectiveness of the present invention or the practicality of the patent. The scope of protection claimed in this application should be determined by the content of its claims, and the specific embodiments described in the specification can be used to interpret the content of the claims.

Claims

1. A three-dimensional inversion method for airborne magnetics based on unstructured node parameterization, characterized by the following steps: include: S1: Define the magnetization amplitude parameter on the nodes of the unstructured tetrahedral mesh, and interpolate it within the tetrahedral element using a linear shape function to obtain a piecewise linear continuous magnetization model. S2: Based on the piecewise linear continuous magnetization model, the analytical response of the linear magnetization of the unit is equivalently recombined to obtain the analytical expression of the forward response consistent with the nodal parameterization. S3: Based on the forward response analytical formula, construct the sensitivity matrix by accumulating the local contributions of the units to the global node columns; S4: Based on the sensitivity matrix, construct an inversion objective function that includes a data fitting term and a model constraint term based on node parameterization; S5: Minimize the norm of the node parameterized model and the norm of the spatial gradient change, solve the constraint terms of the node parameterized model, and thus update the inversion objective function; S6: The Gauss-Newton method is used to iteratively linearize the inversion objective function, and an augmented sparse least squares subproblem is constructed by combining the upper and lower bound constraints of the model parameters with variable transformation. The LSQR algorithm is used to solve the augmented sparse least squares subproblem to obtain the distribution model of underground magnetization amplitude.

2. The three-dimensional inversion method of airborne magnetics based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, The piecewise linear continuous magnetization model is shown below: in, For any spatial location within it The continuous magnetization amplitude at that point, For unit numbering, Numbering of local nodes in the unit. It is a unit The The magnetization amplitude at each local node Representation unit Middle and the first Local scalar linear form functions associated with local nodes; , , , For coefficients; Among them, coefficient , , , As shown below: In the formula, Tetrahedral unit Volume; , , and Tetrahedral unit The geometric coordinates of the four vertices; , , , Number the 4 vertices Positive cyclical replacement.

3. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, The steps for equivalent reconstruction of the analytical response to linear magnetization of a single element include: observation point Set the origin of the local coordinate system, and any source point in this local coordinate system is represented as... ; Construct the magnetization amplitude function, i.e.: in, , Units The constant term and the linear term of the interpolation function; Based on the magnetization amplitude function, a forward-modeling analytical expression consistent with the nodal parameterization is constructed, namely: in, For nodes The magnetization amplitude parameter is defined at this location. As a constant, Permeability in a vacuum environment , For observation point With polyhedra Any source point The distance between them For the first polyhedron individual aspects Outward unit normal vector Observation point to noodles vertical distance, For the number of faces of a polyhedron, Represent a face The total number of edges on, Representing an edge Outward unit normal vector on, This is the gradient operator.

4. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, Sensitivity matrix In the middle, the first Line 1 The column element is the first The observation data for the first The first-order partial derivatives of the node parameters, i.e. in, For nodes When taking the unit amplitude and setting all other nodal parameters to zero, at the observation point The unit magnetic flux density generated at that location; When the observed data is total magnetic anomaly data, the elements of the sensitivity matrix are represented as follows: in For the first The unit vector of the geomagnetic field direction at each observation point.

5. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, The inversion objective function, which includes data fitting terms and node-based parameterization model constraints, is shown below: in, This is the data fit difference term. This represents the norm constraint term of the regularization model based on node parameterization. The regularization factor is used to adjust the weight of the regularization term in the objective function; This represents the magnitude vector of magnetization intensity at the grid nodes.

6. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, The data fit difference term is shown below: in for Diagonal data weight matrix; 1st dimension OK elements Indicates the first Estimate the reciprocal of the noise standard deviation at each observation point; For predicting data vectors; For observation data vectors; The size of the observed data; This is a sensitivity matrix based on node parameterization.

7. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, The steps to solve for the constraint terms of a node-parameterized model include: The model norm constraint term, consisting of the model minimization constraint term and the model spatial gradient minimization constraint term, is expressed as a quadratic stable functional based on the L2 norm, namely: in, For the minimum constraint term of the model based on the L2 norm, The smoothing constraint term for the L2 norm-based node parameterization model is used to penalize the spatial gradient of the magnetization model; Order No. The gradient coefficients of the shape functions corresponding to each local node are: , and The unit volume is Construct the model smoothing constraint matrix as follows: in For unit numbering, Representative Unit The Middle The global number of each node; Calculate the spatial gradient based on the model smoothness constraint matrix. As shown below: in, For the node-parameterized model smoothing constraint matrix, For gradient operators; Based on spatial gradient The model norm constraint terms are calculated.

8. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 1, characterized in that, Before iteratively linearizing the inverted objective function using the Gauss-Newton method, a depth-weighted matrix is ​​used. With horizontal weighted matrix By weighting the node parameterization model, we obtain: in, For nodes depth, The average height of the observation point. The vertical spatial range of the underground inversion model; The value is the vertical magnetization plane observation data or interpolation data corresponding to the underground grid node; Known reference model Introducing regularization constraints, we get: in, , , and These are the weight coefficients of the model's minimum constraint term and the model's smoothness constraint term, respectively.

9. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 8, characterized in that, The augmented sparse least squares subproblem is shown below. in, For the augmented coefficient matrix, To augment the data vector.

10. The three-dimensional airborne magnetic inversion method based on unstructured tetrahedral mesh node parameterization according to claim 9, characterized in that, Before using the LSQR algorithm to solve the augmented sparse least squares subproblem, the augmented sparse least squares subproblem is linearized. The linearized augmented sparse least squares subproblem is shown below: in, , , , These are the upper and lower limits of the model, respectively. This is a diagonal matrix representing the first derivative of the parameter transformation, with diagonal elements of... .

Images