A method and system for 3D spherical shell MT conductivity inversion

Abstract

This invention discloses a method for inverting the three-dimensional magnetotelluric conductivity of a spherical shell, comprising: constructing a three-dimensional finite element mesh of the spherical shell, establishing a mapping relationship between unconstrained parameter vectors and conductivity, and obtaining an initial conductivity model; constructing the frequency domain Maxwell equations based on the model, performing frequency domain forward modeling on each excitation source and frequency after finite element discretization, and solving for the electric field coefficient vector; extracting multi-source tangential electromagnetic field components by interpolation at the observation station locations and through coordinate transformation and Faraday's law derivation; constructing a multi-source impedance tensor using the least squares formula, and defining a total objective function including weighted data fitting terms and regularization terms; solving for the adjoint field vector through the adjoint equation, and calculating the partial derivative of the objective function with respect to the parameter vector; determining the search direction and iteration step size based on the partial derivative, and updating the parameter vector; and outputting the final conductivity model of the spherical shell region after iteration until convergence. This method breaks through the limitations of traditional planar inversion, is suitable for global-scale applications, and improves inversion efficiency and accuracy.

Description

Technical Field

[0001] This invention relates to the field of geophysical inversion technology, particularly to a three-dimensional magnetotelluric (MT) inversion method. Specifically, it utilizes multi-source current excitation and multi-station impedance observation data to invert the conductivity distribution of a spherical crust region through finite element forward modeling and the adjoint gradient method. Background Technology

[0002] Magnetotelluric (MT) is a passive electromagnetic detection method. Researchers infer the electrical conductivity structure of subsurface media by recording changes in the Earth's natural electromagnetic field over time. In other words, we do not actively emit signals, but rather utilize the Earth's own electromagnetic disturbances to obtain subsurface information.

[0003] In practical applications, the application of MT (Mean Transformer) is very broad. It plays an important role in everything from basic science to resource exploration. For example, in global-scale studies, researchers use MT data to infer mantle temperature distribution, water content, and the presence of partial melts; in engineering geophysics, MT is often used to detect mineral resources and geothermal systems; in structural geology studies, it can reveal plate boundaries, fault zones, and deep structural features; and in oil and gas exploration, conductivity structures are often used to identify electrical anomalies in deep reservoirs.

[0004] Currently widely used MT inversion software, such as OCCAM, Rhoplus, and ModEM, are largely similar in their basic principles. Most systems assume the incident electromagnetic field is a plane wave and establish electromagnetic field equations in Cartesian coordinates. The computational object is typically a horizontally layered model or a two- or three-dimensional conductivity structure at a regional scale. To simplify the problem, these methods often employ the assumption of single-source, single-polarization excitation and solve for model parameters using Tikhonov regularization combined with damped least squares iteration. At regional scales (typically less than approximately 500 km), these methods have been extensively validated in practice and have formed a relatively mature technical system.

[0005] As the research scope expands to a global scale, the aforementioned methods gradually reveal significant limitations. The first is the issue of coordinate systems and geometric description. Traditional MT inversions are typically based on the plane wave assumption and modeled in a Cartesian coordinate system. However, on a global scale, the influence of the Earth's curvature can no longer be ignored. Observation stations are distributed across the sphere, and subsurface anomalies also exist within the spherical shell structure. If a planar approximation is still used, it becomes difficult to accurately describe the true geometric relationships. Therefore, these methods often produce significant errors when processing data spanning large geographical areas.

[0006] Another problem stems from the electromagnetic excitation itself. The global electromagnetic field is not a simple single-source, single-polarization structure, but a complex system formed by the superposition of multiple sources, such as magnetospheric currents, ionospheric currents, and tidal magnetic fields. In this context, the applicability of the traditional formula (Z=E / H, where Z is the scalar impedance, E is the scalar value of the horizontal tangential electric field observed at the surface, and H is the scalar value of the horizontal tangential magnetic field orthogonal to the electric field) is limited because it implicitly assumes a single excitation source. For multi-source, multi-polarization electromagnetic fields, there is currently a lack of a unified and clear method for defining the impedance tensor and for its stable calculation. This makes it difficult for existing techniques to fully utilize the information contained in the observational data during actual inversion processes.

[0007] The same problems exist in calculating parameter sensitivity. If inversion is performed in a spherical coordinate system, the observation stations often do not correspond to node positions in the finite element or finite difference mesh, thus requiring three-dimensional spatial interpolation to obtain the electric and magnetic field values. Furthermore, the derivative of the multi-source impedance tensor with respect to the conductivity parameter involves a complex chain of derivatives, making the calculations cumbersome. If the finite difference method is used to estimate sensitivity by perturbing each perturbation model parameter individually, a large amount of forward modeling calculations are often required. In global-scale models, this computational load is generally unacceptable.

[0008] Furthermore, large-scale parameter optimization also faces numerical stability issues. Subsurface electrical conductivity typically spans several orders of magnitude, roughly ranging from 10... -6 S / m to 10 3 S / m. Such a large numerical range can lead to significant imbalances in the optimization problem. If constraints or regularization strategies are not handled properly, the inversion process may encounter convergence difficulties or even numerical instability.

[0009] From a software implementation perspective, existing spherical coordinate MT forward modeling programs also have certain limitations. Some codes only provide electromagnetic field forward modeling functions and do not form a complete inversion system; the sensitivity of multi-source impedance tensors is usually not systematically derived and implemented; at the same time, publicly available, clearly structured, and easily verifiable global-scale MT inversion software is still relatively limited. Summary of the Invention

[0010] The purpose of this application is to provide a method and system for inverting the three-dimensional magnetotelluric conductivity of a spherical shell, which improves the inversion efficiency and accuracy.

[0011] The technical solution provided in this application is as follows:

[0012] In a first aspect, this application provides a method for inverting the three-dimensional magnetotelluric conductivity of a spherical shell, including:

[0013] S1. Construct a three-dimensional finite element mesh for the spherical shell, set the initial conductivity value of each element in the three-dimensional finite element mesh for the spherical shell, establish the mapping relationship between the unconstrained parameter vector and the conductivity, and obtain the initial conductivity model of the spherical shell region.

[0014] S2. Based on the conductivity model, construct the corresponding frequency domain Maxwell equations; based on the frequency domain Maxwell equations, perform finite element discrete frequency domain forward modeling for each excitation source and frequency, and solve for the electric field coefficient vector.

[0015] S3. Based on the three-dimensional finite element mesh of the spherical shell and the electric field coefficient vector, three-dimensional spatial interpolation is performed at the location of the observation station, and after coordinate transformation and Faraday's law derivation, the multi-source tangential electromagnetic field components of each station are extracted.

[0016] S4. Based on the multi-source tangential electromagnetic field components of each station, the multi-source impedance tensor is constructed by calculating the least squares formula.

[0017] S5. Based on the multi-source impedance tensor, define the overall objective function composed of the weighted data fitting term and the regularization term;

[0018] S6. Based on the derivative of the data fitting term with respect to the electric field coefficient vector, construct the adjoint equation and solve for the adjoint field vector; combine the adjoint field vector to calculate the partial derivative of the objective function with respect to the parameter vector.

[0019] S7. Based on the partial derivative of the objective function with respect to the parameter vector, calculate the search direction and determine the iteration step size, and update the parameter vector;

[0020] S8. Repeat steps S2 to S7 for iterative calculations until the convergence condition is met or the maximum number of iterations is reached, and output the final conductivity model of the spherical shell region.

[0021] In one possible implementation, the mapping relationship between the unconstrained parameter vector and conductivity in step S1 is established using Logistic parameterization, expressed as:

[0022] ;

[0023] The derivative is: ;

[0024] in, This represents the conductivity of the r-th grid cell. and These are the lower and upper physical bounds of conductivity, respectively. It is a parameter vector The first in One parameter.

[0025] In one possible implementation, step S2, when performing finite element discrete frequency-domain forward modeling on each excitation source and frequency based on the frequency-domain Maxwell equations, specifically includes: discretizing the frequency-domain Maxwell equations into a complex linear system using finite element methods. ,in For the system matrix corresponding to the frequency, Let be the electric field coefficient vector corresponding to the s-th excitation source. This is the source term vector corresponding to the s-th excitation source;

[0026] The complex linear system is solved using the preconditioned FGMRES iterative solver to obtain the electric field coefficient vector;

[0027] At the same frequency, the system matrix It is only related to the conductivity model, grid and frequency, and is independent of the excitation source. Therefore, for multiple excitation sources at the same frequency, the system matrix and triangular decomposition results are reused, and only the source term vectors corresponding to each excitation source are updated, which greatly reduces the cost of multi-source forward modeling.

[0028] Parallel computing using a message passing interface (MPI) is employed, with each process handling different excitation sources or grid subdomains to achieve distributed accelerated computing for multi-source forward modeling.

[0029] In one possible implementation, the specific process of extracting the multi-source tangential electromagnetic field components of each station by performing three-dimensional spatial interpolation at the station location and deriving the coordinates and Faraday's law in step S3 is as follows: calculate the centroid coordinates of the tetrahedral element where the station is located, calculate the Cartesian component of the electromagnetic field at the station location using the Nédélec edge element basis function, convert it to the tangential electric field component through the spherical coordinate Jacobi matrix, and then derive the tangential magnetic field component from the electric field curl according to Faraday's law.

[0030] In one possible implementation, the least squares formula for constructing the multi-source impedance tensor in step S4 is:

[0031] ;

[0032] in, Let k represent the magnetotelluric impedance tensor of station k. Indicates conjugate transpose; and The tangential electric field matrix and tangential magnetic field matrix are composed of the multi-source tangential electric field and magnetic field components of station k.

[0033] The above and The k station is composed of electromagnetic fields from multiple sources stacked together. matrix, For the number of sources; compared to a single source The formula improves the robustness of impedance estimation.

[0034] In one possible implementation, the overall objective function in step S5 is:

[0035] ;

[0036] ;

[0037] ;

[0038] in, Let be the overall objective function. For the first The regularization parameter for the next iteration; This represents the objective function for fitting the data. Represents the parameter vector The corresponding conductivity model, the predicted impedance real vector obtained through forward modeling and impedance construction; Represents the measured real vector of impedance; Represents a data weighting matrix; This represents the standard function for regularization. This represents the three-dimensional solution domain of the spherical shell. Indicates spatial location The depth weight at each location (decays exponentially with depth). Represents the parameter vector In spatial location The gradient at that point.

[0039] In one possible implementation, step S6 specifically includes:

[0040] S61. Link the derivatives of the data fitting terms with respect to the impedance tensor to the derivatives with respect to the electric field coefficient vector using the chain rule, construct the adjoint equation, and solve it to obtain the adjoint field vectors corresponding to each excitation source.

[0041] S62. Combine the partial derivatives of the adjoint field vector and the mass matrix with respect to the parameter vector, and calculate the partial derivatives of the data fitting term with respect to the parameter vector using the chain rule.

[0042] S63. Calculate the gradient of the regularization term with respect to the parameter vector, and sum the gradient of the data fitting term and the gradient of the regularization term with weights to obtain the partial derivative of the total objective function with respect to the parameter vector.

[0043] The adjoint equation is in the form of: ;

[0044] In the formula, Represents the conjugate transpose of the system matrix. This represents the adjoint field vector corresponding to the s-th excitation source. This represents the electric field coefficient vector corresponding to the s-th excitation source. This is the partial derivative of the data fitting term with respect to the electric field coefficient vector;

[0045] The data fitting term applies to the r-th parameter in the parameter vector. The partial derivatives are calculated using the adjoint-state method, and the formula is:

[0046] ;

[0047] In the formula, The electromagnetic field angular frequency, The permeability of free space, Operator for extracting the real part; Let mr be the partial derivative of the mass matrix with respect to the r-th parameter, and this partial derivative is non-zero only within the mesh cell corresponding to mr.

[0048] In one possible implementation, in step S7, the search direction is calculated using the L-BFGS algorithm, and the iteration step size is determined using the Armijo line search.

[0049] In one possible implementation, the L-BFGS algorithm is configured to store 5 to 10 pairs of historical vectors to reduce memory requirements, and the Armijo line search coefficients are set to... The iteration step size range is ~ The maximum / minimum step size are respectively and To prevent numerical overflow.

[0050] In one possible implementation, during the iteration of S8, the regularization parameter of the regularization term in step s5 is dynamically adjusted based on the decrease in RMS error after the parameter vector update; wherein the dynamic adjustment adopts an exponential decay strategy, and the expression is:

[0051] ;

[0052] in, These are the initial values ​​for the regularization parameters. This is the lower bound of the regularization parameter. The decay rate is the regularization parameter.

[0053] Secondly, this application provides a spherical shell three-dimensional magnetotelluric conductivity inversion system for implementing the above-mentioned method, including an input module, a forward calculation module, an inversion calculation module, an optimization module, and an output module;

[0054] The input module is used to perform the following: constructing a three-dimensional finite element mesh for the spherical shell, setting the initial conductivity value of each element in the three-dimensional finite element mesh for the spherical shell, establishing the mapping relationship between the unconstrained parameter vector and the conductivity, and obtaining the initial conductivity model of the spherical shell region;

[0055] The forward modeling module is used to perform the following: based on the conductivity model, construct the corresponding frequency domain Maxwell equations; based on the frequency domain Maxwell equations, perform finite element discrete frequency domain forward modeling for each excitation source and frequency, and solve for the electric field coefficient vector.

[0056] Based on the three-dimensional finite element mesh of the spherical shell and the electric field coefficient vector, three-dimensional spatial interpolation is performed at the location of the observation station, and after coordinate transformation and Faraday's law derivation, the multi-source tangential electromagnetic field components of each station are extracted.

[0057] The inversion calculation module is used to perform the following: based on the multi-source tangential electromagnetic field components of each station, the least squares formula is used to calculate and construct the multi-source impedance tensor.

[0058] Based on the multi-source impedance tensor, a total objective function composed of a weighted data fitting term and a regularization term is defined;

[0059] Based on the derivative of the data fitting term with respect to the electric field coefficient vector, the adjoint equation is constructed and the adjoint field vector is obtained by solving it; combined with the adjoint field vector, the partial derivative of the objective function with respect to the parameter vector is calculated.

[0060] The optimization module is used to perform the following: calculate the search direction and determine the iteration step size based on the partial derivative of the objective function with respect to the parameter vector, and update the parameter vector; perform iterative calculations through the forward calculation module and the inverse calculation module until the convergence condition is met or the maximum number of iterations is reached;

[0061] The output module is used to output the final conductivity model of the spherical shell region.

[0062] In one possible implementation, the system further includes a parallel module for implementing parallel computing.

[0063] In the above scheme, the system reads basic data through the input module, performs frequency domain forward modeling and electromagnetic field interpolation through the forward modeling module, constructs the impedance tensor and calculates the parameter gradient through the inversion modeling module, implements parameter updates and regularization adjustments through the optimization module, determines convergence and outputs the conductivity model through the output module, and implements MPI parallelism, global gradient reduction and load balancing through the parallel module.

[0064] Thirdly, this application provides an electronic device, including: a memory and a processor;

[0065] The memory is used to store computer programs;

[0066] The processor is used to invoke the computer program to execute the method described above.

[0067] Fourthly, this application provides a computer-readable storage medium storing a computer program that, when executed on an electronic device, causes the electronic device to perform the method described above.

[0068] Fifthly, this application provides a computer program product, including a computer program that, when run on an electronic device, causes the electronic device to perform the method described above.

[0069] The specific implementation methods of the second to fifth aspects of this application can refer to the implementation methods of the first aspect, and will not be elaborated here.

[0070] The method and system for inverting the three-dimensional magnetotelluric conductivity of a spherical shell provided in this application have the following advantages:

[0071] Global scale adaptability: For the first time, a magnetotelluric inversion system with a three-dimensional finite element mesh for the spherical shell is constructed, breaking through the limitations of the local assumptions of traditional planar inversion. It can be directly applied to global scale subsurface conductivity imaging, covering the global crust region from the surface to the depths, and adapting to the needs of global geophysical exploration.

[0072] Multi-source impedance robustness: By constructing a multi-source impedance tensor using multi-source tangential electromagnetic field components, the noise resistance and reliability of the impedance tensor are greatly improved compared to traditional single-source impedance estimation, solving the technical problems of single-source observation being susceptible to noise interference and impedance estimation distortion.

[0073] The efficiency of the adjoint gradient: By constructing a mathematically consistent adjoint gradient system, in each round of inversion iteration, only one multi-source multi-frequency forward modeling + one multi-source adjoint equation solution is required to obtain the full parameter gradient. Compared with the finite difference method (one forward modeling per parameter), the efficiency is improved by several orders of magnitude, making it possible to perform large-scale inversion of hundreds of thousands of parameters under the global spherical shell grid.

[0074] Multi-source forward modeling acceleration: Multiple excitation sources at the same frequency reuse the system matrix and triangular decomposition results, only updating the source term vector. Combined with parallel computing, the computational cost of multi-source forward modeling is greatly reduced, realizing distributed parallel acceleration at the source level and frequency level.

[0075] Inversion stability and accuracy: An overall objective function combining an unconstrained parameter-conductivity logistic mapping, a weighted data fitting term, and a Laplace regularization term is adopted, along with L-BFGS optimization and Armijo line search, to ensure stable convergence of the inversion iteration and improve the recovery accuracy and spatial resolution of the conductivity model.

[0076] Technical scalability: For the first time, the adjoint method is extended to multi-source tensor observations (magnetic impedance) in spherical coordinates. The complete chain derivation of parameter gradients from tensor observations is realized through tools such as the Kronecker product, providing a general technical framework for subsequent multi-source tensor electromagnetic inversion. Attached Figure Description

[0077] Figure 1 This application presents a flowchart of a method for inverting the three-dimensional magnetotelluric conductivity of a spherical shell according to an embodiment.

[0078] Figure 2 This application verifies the truth diagram of the chessboard model in Case 1, corresponding to a spherical shell slice at a depth of 10km.

[0079] Figure 3 The diagram shows the inversion evolution process of the chessboard model in Case 1 of this application, corresponding to a spherical shell slice with a depth of 100km. (a), (b), (c), (d), (e), and (f) represent the inversion results for iterations of 1, 10, 30, 50, 100, and 200 steps, respectively.

[0080] Figure 4 The convergence curve of the chessboard model in Case 1 of this application is shown in (a), (b), (c), (d), (e), and (f), respectively. The curves show the trends of data fitting term, regularization term, regularization parameter, root mean square error, gradient norm, and step size as a function of the number of iterations.

[0081] Figure 5 This application verifies the true value of the local block anomaly model in Case 2, which corresponds to a spherical shell slice view at a depth of 10km.

[0082] Figure 6 The following diagram illustrates the inversion evolution process of the local block anomaly model in Case 2 of this application, where (a), (b), (c), (d), (e), and (f) represent the inversion results corresponding to iteration 1, iteration 5, iteration 10, iteration 30, iteration 50, and iteration 100, respectively.

[0083] Figure 7 The convergence curve of the local block anomaly model in Case 2 of this application is shown in (a), (b), (c), (d), (e), and (f), respectively. The curves show the trends of data fitting term, regularization term, regularization parameter, root mean square error, gradient norm, and step size with the number of iterations.

[0084] Figure 8 This application presents a schematic diagram of the module architecture of a spherical shell three-dimensional magnetotelluric conductivity inversion system according to an embodiment. Detailed Implementation

[0085] To enable those skilled in the art to better understand the present application, the technical solution of the present application will be further described in detail below with reference to the embodiments and accompanying drawings.

[0086] Example 1:

[0087] like Figure 1 As shown, this embodiment provides a method for inverting the three-dimensional magnetotelluric conductivity of a spherical shell. The overall technical framework is to construct an integrated three-dimensional inversion system for spherical shells, encompassing forward modeling, sensitivity assessment, and optimization. This system overcomes the limitations of the Cartesian coordinate system in traditional planar magnetotelluric (MT) inversion and adapts to global-scale spherical shell geometric features and multi-source electromagnetic excitation scenarios. The method consists of four core steps, each interconnected: from basic mesh construction and parameter initialization, to multi-source electromagnetic field forward modeling and impedance tensor construction, to efficient parameter sensitivity calculation, and finally, to iterative updates and convergence of the conductivity model through optimization algorithms. The complete process covers the entire chain from data input to inversion result output.

[0088] To enable those skilled in the art to better understand the present application, the following detailed description of the specific implementation method of the present application is provided in conjunction with the technical principles, process design, and practical verification cases of spherical coordinate system magnetotelluric inversion. This implementation method revolves around a four-layer progressive technical framework: spherical shell mesh and parameter initialization, multi-source forward modeling and impedance construction, sensitivity and gradient calculation, L-BFGS optimization and adaptive regularization. At the same time, the effectiveness and accuracy of the method are verified through actual model inversion.

[0089] First layer: Spherical shell mesh and parameter initialization.

[0090] This step lays the foundation for the inversion calculation. The core is to construct a three-dimensional finite element mesh that adapts to the Earth's crust structure, complete the initialization of the conductivity parameters and the parameterized mapping of bounded constraints, and ensure that the initial inversion values ​​conform to physical laws and that the parameter solutions are within a reasonable domain.

[0091] C11. Create a three-dimensional finite element mesh for the spherical shell:

[0092] In some embodiments, the method includes: generating a three-dimensional tetrahedral unstructured mesh of a spherical shell using a three-dimensional finite element mesh generator (such as Gmsh). The mesh covers the Earth's surface (radial radius R = 6371 km) to a specified underground depth (such as 700 km, corresponding to a radial radius R = 5670 km). The outer and inner boundaries of the mesh are both spherical surfaces. The spherical center angle boundary condition is used to adapt to the geometric characteristics of the spherical coordinate system. The number of mesh elements can be adjusted according to the inversion accuracy requirements, and the typical scale can reach hundreds of thousands.

[0093] C12. Coordinate System Setting and Transformation:

[0094] Physical coordinates: spherical coordinates are used. ,in, Radial (from the Earth's center). and These are the radial radii of the inner and outer boundaries of the spherical shell, respectively. The coordinates are relative to the North Pole (based on the North Pole). Longitude.

[0095] Calculated coordinates: To adapt to finite element discrete calculations, spherical coordinates are converted to Cartesian coordinates. The conversion formula is:

[0096] ;

[0097] In the formula, This represents the three components of a Cartesian coordinate system.

[0098] Station coordinate transformation: The original coordinates of the observation stations are spherical coordinates, which need to be uniformly converted into physical coordinates of the finite element mesh to prepare for subsequent electromagnetic field interpolation of the station locations.

[0099] C13. Initial Model and Parametric Mapping of Electrical Conductivity:

[0100] Initial conductivity model: A uniform background model can be used to set a uniform initial conductivity value for all elements of the spherical shell 3D finite element mesh, which serves as the initial benchmark for inversion iteration and ensures the stability of the iteration starting point.

[0101] Logistic parameterization mapping: to satisfy the physical constraints of conductivity ( (Conductivity is non-negative and bounded), establish an unconstrained parameter vector. With conductivity The mapping relationship transforms the constrained conductivity inversion problem into an unconstrained parameter optimization problem. The mapping formula and derivative are as follows:

[0102] ;

[0103] In the formula, This represents the conductivity of the r-th grid cell. and These are the lower and upper physical bounds of conductivity, respectively. parameter vector The r-th parameter in , express dimensional real number field, This represents the total number of model parameters.

[0104] The parametric derivative is:

[0105] .

[0106] The advantages of this mapping are: it achieves the unification of unconstrained parameterization and bounded constraints, automatically satisfies the physical constraints of conductivity (non-negative and bounded), the derivative is smooth, which facilitates gradient calculation, the objective function has better convexity after parameterization, and it is suitable for stable inversion of large-scale parameters (tens of thousands to hundreds of thousands of parameters).

[0107] Second layer: Multi-source forward modeling and impedance construction.

[0108] This step is the core forward modeling component of the inversion process, based on Maxwell's equations in the frequency domain, and utilizes the Sobolev space... Finite element discretization is used to realize the forward modeling of electromagnetic fields under multi-source current excitation. Then, through station electromagnetic field interpolation and multi-source least squares formula, a more robust magnetotelluric impedance tensor is constructed to make full use of the information advantage of multi-source excitation.

[0109] For each frequency and source, the electric field is obtained by solving Maxwell's equations in the frequency domain. Three-dimensional spatial interpolation is performed at the observation station location to extract the tangential electromagnetic field. Multi-source electromagnetic field data are stacked, and the impedance tensor is calculated using the least squares formula. Real number block processing is used to improve numerical stability.

[0110] C21. Governing equations and quasi-static approximation:

[0111] Governing equations:

[0112] To address the low-frequency characteristics of magnetotelluric inversion, neglecting displacement current, the frequency-domain Maxwell's equations of the quasi-static approximation are adopted as the governing equations, the strong form of which is:

[0113] ;

[0114] in, For curl operator, H / m is the permeability of free space. For the frequency domain complex electric field intensity, The imaginary unit, Angular frequency, For frequency; This represents the spatial conductivity distribution (inversion core parameters). Excitation is based on volume current density; and Three-dimensional spatial position vector The function reflects the spatial distribution characteristics of the electrical conductivity of the geological body and the excitation source.

[0115] C22. Weak form derivation of the equation:

[0116] To adapt to finite element numerical solutions, the strong form of the frequency domain Maxwell's equations is transformed into... The variational weak form under the following conditions, for all belonging to trial functions of space ,satisfy:

[0117] ;

[0118] In the formula, For the three-dimensional solution domain of the spherical shell inversion, For trial functions .

[0119] The weak form transforms partial differential equations into integral equations, providing a mathematical foundation for finite element discretization, and simultaneously... Spatial constraints ensure the continuity of the tangential component of the electric field and avoid numerical pseudo-solutions.

[0120] C23 Finite element discretization:

[0121] use Nédélec Type I boundary element basis functions Discretizing the electric field, the basis function is defined only on the grid edges, perfectly fitting the physical property of tangential continuity of the electric field. The discretization formula is as follows:

[0122] ;

[0123] In the formula, The finite element approximation of the frequency domain complex electric field intensity is used to replace the analytical solution. ; This represents the electric field coefficient corresponding to the t-th degree of freedom. Let the t-th Nédélec type-1 basis function satisfy... The requirement for spatial continuity; The total number of degrees of freedom of the spherical shell finite element mesh is equal to the total number of edges of the mesh, corresponding to the dimension of the electric field coefficient vector.

[0124] After discretization, the equations of the complex linear system are obtained:

[0125] ;

[0126] in, This represents the frequency domain complex system matrix (also simply called the system matrix), with dimension 1. ; This represents the electric field coefficient vector, with dimension . , is the core solution objective of finite element forward modeling, and its t-th element is . ; The right-hand vector (source term vector) of the complex linear system equations is derived from... The discrete excitation vector has a dimension of .

[0127] The system matrix , To and The relevant quality matrix, Here is the stiffness matrix, with dimensions equal to... Consistent. The element in the i-th row and t-th column for: ; The element in the i-th row and t-th column for: ;

[0128] The i-th component for: .

[0129] C24. Multi-source frequency domain forward modeling solution:

[0130] To improve the efficiency of forward modeling under multi-source excitation, a pre-conditioned Flexible Generalized Minimal Residual (FGMRES) iterative solver is used to solve the discretized complex linear system. Simultaneously, real-number block processing is employed to enhance numerical stability. Multi-CPU / multi-node collaborative distributed parallel computing technology (MPI parallel computing) enables distributed processing of large-scale computations. The specific solution process is as follows:

[0131] Matrix Assembly and Preconditioning (One-Time): Assembling the Stiffness Matrix from the Mesh and Conductivity Model quality matrix Forming a system matrix By utilizing the Auxiliary Space Maxwell Solver (AMS) of the Modular Finite Element Library (MFEM) to construct precondition operators, iterative convergence is accelerated.

[0132] Multi-source iterative solution (source-by-source): For each excitation source (For different polarizations or periods), read the volume current density from the source file. and cycle Calculate the source term vector corresponding to the excitation source. The linear system is solved using the FGMRES solver, with the convergence condition being: .in, This represents the period of the electromagnetic field corresponding to the s-th excitation source. This represents the electromagnetic field angular frequency corresponding to the s-th excitation source. This represents the residual vector during iterative solution of the s-th excitation source. This represents the convergence tolerance of the FGMRES iterative solution, a mathematical constant (usually taken as...). ).

[0133] In some embodiments, for scenarios with multiple excitation sources at the same frequency, the system matrix is ​​completed in one go during the matrix assembly stage. The triangular decomposition (LU decomposition) is used; in subsequent multi-excitation source solutions, only the corresponding source term vector needs to be updated. The electric field coefficient vector is solved quickly by reusing the completed LU decomposition results through a forward-backward substitution. This eliminates the need to repeatedly perform matrix factorization, significantly reducing the computational cost of multi-source forward modeling.

[0134] The system matrix and LU decomposition results of multiple excitation sources at the same frequency are reused, and only the source term vector is updated, which greatly reduces the computational cost.

[0135] Real-number block processing (to ensure numerical stability): The complex linear system equations are transformed into equivalent real-block system equations, avoiding numerical instability problems in complex numerical computation. The real-block system equations are in the following form:

[0136] ;

[0137] in, Represents the mass matrix The real part, Represents the mass matrix The imaginary part; Let represent the real part of the electric field coefficient vector corresponding to the s-th excitation source. This represents the imaginary part of the electric field coefficient vector corresponding to the s-th excitation source. Let represent the real part of the right-hand vector corresponding to the s-th excitation source. Let represent the imaginary part of the right-hand vector corresponding to the s-th excitation source.

[0138] C25. Tangential electromagnetic field interpolation at the station:

[0139] Not all observation stations are located at the nodes of the finite element mesh. Therefore, it is necessary to extract the electromagnetic field at the station location through three-dimensional spatial interpolation, and then obtain the tangential electric and magnetic field components of the station using coordinate transformation and Faraday's law. The specific process is as follows:

[0140] Definition of tangential coordinate system: at the station location The tangential basis vector in spherical coordinates is (North-South) and (East-west direction).

[0141] The Jacobian matrix corresponding to the Cartesian coordinate system is:

[0142] ;

[0143] in, , , These are the radial unit basis vector, the north-south tangential unit basis vector, and the east-west tangential unit basis vector in spherical coordinates, respectively. The covariance of the spherical coordinate system. Longitude in spherical coordinates. , , These are the unit basis vectors in the x, y, and z directions of the Cartesian coordinate system, respectively.

[0144] Electric field interpolation: Calculate the barycentric coordinates of the tetrahedral element containing the station, and use Nédélec edge element basis functions to calculate the electric field vector of the station's location in Cartesian coordinates. The formula is:

[0145] ;

[0146] in, Let k be the spatial position vector of the kth observation station. The index for the observation stations is a mathematical traversal variable that iterates through all surface observation stations, with a range of [missing information]. , Total number of stations; This represents the finite element approximation of the complex electric field in the frequency domain at the location of the k-th observation station. Indicates the basis function of the t-th Nédélec type I edge element in The value at that location.

[0147] Then, using the spherical coordinate Jacobian matrix (rotation matrix) The electric field vector in the Cartesian coordinate system Tangential electric field components converted to spherical coordinates . ,in and spherical coordinates Perpendicular to the radial direction The components are the north-south tangential electric field component and the east-west tangential electric field component, respectively.

[0148] Magnetic field interpolation: The tangential magnetic field component is derived from the curl of the electric field according to Faraday's law, and the magnetic induction intensity is... The calculation formula is:

[0149] ;

[0150] In the formula, This represents the magnetic flux density vector.

[0151] The approximate solution for the magnetic flux density at the station location under the finite element framework is:

[0152] ;

[0153] In the formula, This represents the finite element approximation of the magnetic induction intensity at the location of station k. express The curl of the value is taken at the k-th station.

[0154] Then, the magnetic field strength in the Cartesian coordinate system is converted into a magnetic field strength vector in the spherical coordinate system using the spherical coordinate Jacobian matrix, and then... Converted to tangential magnetic field components .

[0155] C26. Construction of the Multi-Source Impedance Tensor

[0156] This application overcomes the limitations of the traditional single-source impedance formula Z=E / H under multi-source excitation. It employs a multi-source least squares formula to construct the impedance tensor, fully utilizing information from multiple independent excitations to improve the robustness of impedance estimation and automatically handle polarization ambiguity issues. Specifically, it includes:

[0157] Multi-source electromagnetic field matrix stacking: for site ,Will The tangential electric and magnetic field components under each excitation source are stacked into a matrix, in the form of:

[0158] ;

[0159] ;

[0160] In the formula, This indicates that station number k is located at... Tangential electric field matrix under single-source excitation This indicates that station k is located at position 1~ Tangential electric field components under individual source excitation; This indicates that station number k is located at... Tangential magnetic field matrix under individual source excitation This indicates that station k is located at position 1~ Tangential magnetic field components under individual source excitation; Represents the complex field.

[0161] The formula for calculating the least-squares impedance tensor of multiple sources is as follows:

[0162] ;

[0163] in, Let k represent the magnetotelluric impedance tensor of station k. This indicates the conjugate transpose.

[0164] This formula is based on the least squares principle, making full use of the information provided by multiple independent excitation sources, improving the robustness of impedance estimation, and automatically handling polarization confusion problems.

[0165] Impedance observation and objective function:

[0166] Real vectorization of impedance tensor: The impedance tensor, when flattened by rows or columns, is a 4-dimensional complex vector, as shown in the formula: , For the 2×2 impedance tensor of station k The four complex components are then stacked according to their real and imaginary parts to form an 8-dimensional real vector, in the form of:

[0167] ;

[0168] In the formula, The real-vectorized vector representing the impedance of station k. The real part extraction operator is used to represent the real part extraction operator. The real part extraction operator is used to represent the real part extraction operator. It represents an 8-dimensional real space. The real-vectorized impedance provides the numerical basis for subsequent objective function construction and gradient calculation.

[0169] Third layer: Sensitivity and gradient calculation

[0170] Parameter sensitivity is the core of inversion optimization. This step uses the adjoint gradient method to calculate the sensitivity of the conductivity parameter, overcoming the limitation of large computational cost in the traditional finite difference method. By deriving the analytical Jacobian matrix of the impedance tensor with respect to the electric and magnetic fields, and combining it with the chain rule, efficient calculation of the parameter gradient from impedance observation is achieved. The gradient calculation efficiency is improved compared to the finite difference method. This multiplier makes it possible to retrieve large parameters on a global scale.

[0171] C31. Construction of the overall objective function:

[0172] Construct a data objective function and use weighted average. The norm form measures the deviation between the predicted impedance and the measured impedance, and the formula is:

[0173] ;

[0174] In the formula, This represents the objective function for fitting the data. Represents an unconstrained parameter vector The corresponding conductivity model, the predicted impedance real vector obtained through forward modeling and impedance construction; Represents the measured real vector of impedance; This represents a weighted matrix of data.

[0175] Constructing the regularization term: In some embodiments, depth-weighted regularization in the Günther (2006) format is used to constrain the spatial gradient of the model parameters, ensure model smoothness, and suppress spurious anomalies. The formula is as follows:

[0176] ;

[0177] In the formula, This represents the standard function for regularization. This represents the three-dimensional solution domain of the spherical shell. Indicates spatial location Depth weight at the location, Represents the parameter vector In spatial location The gradient at that point.

[0178] The discretized form is: ;

[0179] in, The model regularization weighting matrix consists of the grid Laplacian operator and depth weights.

[0180] Overall objective function:

[0181] ;

[0182] in, Let be the overall objective function. For the first The regularization parameter for the next iteration.

[0183] C32. Derivation of the impedance Jacobian matrix:

[0184] The impedance Jacobian matrix reflects the partial derivatives of the impedance tensor with respect to the electric and magnetic fields, and is the core mathematical foundation of the adjoint gradient method. For multi-source least squares impedance, let the magnetic field Gram matrix of the k-th station be... ,but Derive its Jacobian matrices with respect to the electric and magnetic fields respectively:

[0185] The Jacobian matrix of the electric field (in vectorized form) is:

[0186] ;

[0187] in, Represents vectorization operators; It is the Kronecker product. It is a 2-order identity matrix.

[0188] For the Jacobian matrix of the magnetic field: the derivation of the derivative of the matrix inverse is involved, and the partial derivative is calculated by combining the chain rule, which provides a basis for the subsequent construction of the adjoint source.

[0189] C33. Construction and solution of the adjoint equation:

[0190] To calculate the gradient of the data fitting term with respect to the conductivity parameter, we first use the chain rule to link the derivative of the data fitting term with respect to the impedance tensor to the electric field coefficient vector. The derivative of is used to construct the adjoint equation and solve for the adjoint field vector. The adjoint equation is of the form:

[0191] ;

[0192] In the formula, Represents the conjugate transpose of the system matrix. This represents the adjoint field vector corresponding to the s-th excitation source. This represents the electric field coefficient vector corresponding to the s-th excitation source. This is the right-hand vector of the adjoint equation (residual term), which is the partial derivative of the data fitting term with respect to the electric field coefficient vector.

[0193] The gradient of the regularization term with respect to the parameters can be directly calculated from the gradient in the model space without solving the adjoint equation; the gradient of the final objective function is the weighted sum of the gradient of the data fitting term and the gradient of the regularization term.

[0194] C34. Parameter sensitivity (gradient) calculation:

[0195] By combining the solved adjoint field vector and the partial derivatives of the mass matrix with respect to the parameters, the chain rule is used to complete the data fitting term with respect to the parameter vector. The gradient of the objective function is calculated, and then the gradient of the regularization term is superimposed to obtain the complete gradient of the objective function with respect to the parameter vector.

[0196] Basic derivation of the chain rule: For the r-th model parameter controlling conductivity That is, parameter vector For the r-th element, the partial derivative of the data fitting term with respect to this parameter can be obtained using the chain rule:

[0197] ;

[0198] In the formula, This represents the partial derivative of the data fitting term with respect to the r-th model parameter. Let represent the partial derivative of the product of the system matrix and the electric field coefficient vector of the s-th source with respect to the r-th parameter.

[0199] Formula for calculating gradient using the adjoint-state method:

[0200] Based on the finite element discretization form of Maxwell's equations in the frequency domain, the above integral formula is simplified into an efficient computational form using the adjoint-state method:

[0201] ;

[0202] in, For the mass matrix The partial derivative with respect to the r-th parameter is only related to the partial derivative with respect to the r-th parameter. The corresponding grid cell is non-zero, while the other cells are zero, which greatly reduces the computational cost and complexity of gradient calculation.

[0203] The standard adjoint method has been used for scalar observations. This invention extends it to multi-source tensor observations (magnetic impedance tensor) in spherical coordinates for the first time. The key is to derive the adjoint source from the 2×2 complex impedance tensor and realize the complete chain derivative derivation from the impedance tensor to the electric field and then to the conductivity parameter through advanced linear algebra tools such as the Kronecker product, thus constructing a mathematically self-consistent adjoint gradient system.

[0204] This method requires only one forward modeling and one solution of the adjoint equation per iteration to obtain the gradient of all parameters. Compared with the finite difference method (one forward modeling per parameter), it improves computational efficiency and makes it possible to perform large-scale inversion of hundreds of thousands of parameters under the global spherical shell grid, realizing efficient chain-like calculation from impedance observation to parameter gradient.

[0205] Fourth layer: L-BFGS optimization and adaptive regularization:

[0206] This step is the optimization iteration stage of the inversion. The finite-memory BFGS (L-BFGS) quasi-Newton algorithm is used to solve the inversion optimization problem, realizing the optimization update of the parameter vector without constraints. Combined with the adaptive regularization decay strategy, the regularization parameter is dynamically adjusted according to the RMS error to achieve "early stability and later refinement" of the model. Finally, the iteration converges to the optimal conductivity model.

[0207] The L-BFGS algorithm does not require storing the complete Hessian matrix; it only stores the parameter update vectors and gradient difference vectors from the most recent k pairs (k is typically 5-10) of iterations. It iteratively constructs an approximate Hessian inverse matrix, significantly reducing memory usage and computational complexity for large-scale model inversion while ensuring superlinear convergence. It is suitable for 3D magnetotelluric conductivity inversion scenarios in a spherical coordinate system. The core steps are as follows:

[0208] Calculate the search direction: Based on the approximate Hessian inverse matrix and the current gradient, calculate the search direction for the nth iteration using the following formula:

[0209] ;

[0210] In the formula, Let be the search direction vector for the nth iteration. Let be the approximate Hessian inverse matrix of the nth iteration. Let be the gradient vector of the overall objective function in the nth iteration. , This represents the unconstrained parameter vector for the nth iteration.

[0211] Line search determines the step size: Along the search direction, the optimal step size is found using the Armijo conditions to ensure that the objective function monotonically decreases. The Armijo conditions are:

[0212] ;

[0213] In the formula, This represents the step size of the nth iteration. For the Armijo line search coefficients, mathematical constants (which can be taken as...) Step size range is limited to ~ To prevent numerical overflow.

[0214] Parameter vector update: Based on the search direction and step size, update the parameter vector for the (n+1)th iteration, using the following formula:

[0215] ;

[0216] In the formula, This represents the unconstrained parameter vector for the (n+1)th iteration.

[0217] Approximate Hessian inverse matrix update: update the parameters using the parameters from the current iteration. and gradient difference Update the approximate Hessian inverse matrix This prepares for the next iteration.

[0218] Adaptive regularization decay strategy:

[0219] Regularization parameters Directly affecting the smoothness and detail resolution of the model, an exponential decay strategy is adopted to dynamically adjust the parameters based on the rate of decrease in RMS error, eliminating the need for manual parameter tuning and achieving automated optimization. The formula is as follows:

[0220] ;

[0221] in, These are the initial values ​​for the regularization parameters. This is the lower bound of the regularization parameter. The decay rate is a regularization parameter, typically ranging from 0.5 to 0.9.

[0222] The core logic of this attenuation strategy is as follows:

[0223] Early iterations: Use strong regularization (larger) (value), suppressing drastic fluctuations in the model, ensuring inversion stability, and rapidly reducing RMS error;

[0224] Mid-term iteration: As the RMS error decreases, the regularization parameter gradually decays, and model details begin to emerge;

[0225] Later iterations: The regularization parameter is reduced to near the lower bound, the model smoothness constraint is weakened, the local details of the model are refined, and the inversion resolution is improved.

[0226] Convergence criteria:

[0227] The iterative process stops when any of the following conditions are met to ensure the accuracy and efficiency of the inversion results. The convergence condition is:

[0228] (1) Gradient norm convergence: The infinite norm of the gradient is sufficiently small. ,in, Represents the infinite norm operator, This represents the gradient convergence tolerance, which is set according to the inversion requirements.

[0229] (2) RMS error convergence: The RMS error is reduced to the target value. ,in This represents the total dimension of the observed data. Indicates the target RMS error;

[0230] (3) Reach the maximum number of iterations: To prevent iteration from getting stuck in an infinite loop, set the maximum number of iterations (e.g., 200 times).

[0231] The converged unconstrained parameter vector The conductivity values ​​of each element in the spherical shell grid are restored using the Logistic parametric mapping formula, resulting in the final conductivity model of the spherical shell region in global spherical coordinates, thus completing the entire inversion process.

[0232] It should be understood that the step numbers S1 to S8 are only used to distinguish and facilitate the expression of different steps / modules, and do not necessarily constitute a restriction on the execution order between the steps.

[0233] The spherical shell three-dimensional magnetotelluric conductivity inversion method provided in the above embodiments of this application calculates the electromagnetic field on a spherical shell three-dimensional finite element mesh based on multi-source volume current excitation using H(curl) finite element discrete frequency domain Maxwell's equations; performs three-dimensional spatial interpolation at the observation station locations, and constructs the impedance tensor using the multi-source least squares formula; uses impedance observation data as constraints, and effectively calculates parameter sensitivity using the impedance Jacobian matrix and the adjoint gradient method; finally, completes the spherical shell conductivity inversion using the adaptive regularized attenuation L-BFGS optimization method. This method establishes a complete multi-source impedance inversion system in spherical coordinates.

[0234] Compared to existing technologies, this invention establishes a complete theoretical system for spherical coordinate MT inversion: defining a multi-source impedance tensor within a spherical shell finite element framework and effectively calculating its sensitivity to conductivity. Breaking through the geographical limitations of traditional planar MT inversion, the multi-source impedance tensor provides richer information than single-source observations, improving resolution and stability. Automatic coupling of unconstrained optimization and physical constraints is achieved through Logistic parameterization; a complete chain rule from impedance observation to parameter gradient is derived; and an adaptive regularization attenuation strategy is designed to automatically adjust according to iterative progress. This integrates forward modeling, sensitivity assessment, and optimization into a scalable inversion system, supporting efficient inversion of large-scale data from multiple sources and stations simultaneously, supporting global-scale applications, breaking the geographical limitations of traditional planar MT inversion, and opening new directions for global mantle conductivity imaging. It supports applications such as global mantle temperature structure, water content distribution, and melt detection.

[0235] Actual model inversion verification

[0236] To verify the effectiveness, accuracy, and robustness of this method, two typical models, a chessboard model and a local block anomaly model, were designed for inversion verification, covering two geological scenarios: periodic anomalies and local concentrated anomalies. The inversion effect was evaluated from dimensions such as data fitting accuracy, model feature recovery, and spatial positioning accuracy.

[0237] Verification Case 1: Impedance Inversion of the Chessboard Model.

[0238] Model and data preparation:

[0239] Grid and Truth

[0240] Spherical grid: Earth's surface up to a depth of ~700km It has a total of 131,842 tetrahedral elements and approximately 100,000 degrees of freedom.

[0241] Initial conductivity: Initial uniform background S / m, the true conductivity value is a checkerboard model with a period of 150km (a periodic conductivity distribution model (high conductivity and low conductivity regions are arranged alternately like a checkerboard), such as Figure 2 The image shown is a slice view of the spherical shell at a depth of 10 km. scope S / m.

[0242] Excitation sources and observations: 3 different polarization excitation sources, periodic Second( Hz, (rad / s); 1200 stations are evenly distributed on the surface. The observation data are complex values ​​with 4 impedance components, with a total dimension of 9600 (8 real and imaginary components × 1200 stations), and 1% Gaussian noise is added to simulate actual observation errors.

[0243] Inversion parameter settings:

[0244] Parameterization: One parameter per unit Logistic mapping is used. The initial values ​​are the parameters corresponding to the uniform background.

[0245] Regularization: Günther (2006) depth-weighted, initial value , , .

[0246] L-BFGS memory Maximum iterations 200, gradient tolerance Target RMS error =0.01.

[0247] Computing resources: 32 MPI processes.

[0248] Inversion results and analysis:

[0249] Model feature recovery: Figure 3 The diagram shows the corresponding inversion evolution process. This process clearly demonstrates the gradual identification and evolution from an initial uniform background to the final chessboard model. The large-scale framework of the chessboard emerges in the mid-stage of the iterations, while the structure becomes clearer in the later stages, with local details highly consistent with the ground truth. This evolution process clearly demonstrates the inversion's ability to gradually identify and recover periodic structures.

[0250] Iterative Convergence Characteristics: The iterative process (summary) is shown in Table 1. As shown in Table 1, in the early stage (iterations 1 and 20): RMS rapidly decreased from 38.7 to 15.96, the objective function decreased by 93%, and the search direction became clear. In the middle stage (iterations 21 and 100): RMS decreased slowly, regularization gradually decayed, and chessboard details became apparent. In the later stage (iterations 101 and 200): RMS tended to 0.05~0.1, eventually converging to an accuracy of <1%. The RMS error monotonically decreased from the initial 38.69 to 0.052, and the data fitting accuracy improved by 99.9%; the objective function... The objective function shows a stable decreasing trend as it decreases from 1.94 × 10⁸ to 2.06 × 10⁴, as shown below. Figure 4 As shown, the convergence curve exhibits a stable characteristic of monotonically decreasing. The regularization term makes a significant contribution in the early stages, stabilizes in the later stages, and shows no divergence in the inversion process.

[0251] Table 1. Iterative Process (Summary) for Validation Case 1

[0252] ;

[0253] Accuracy and positioning: The true conductivity is about 0.1 S / m in the high conductivity region and about 0.001 S / m in the low conductivity region. The conductivity range of the inversion result is close to the true value, and the contrast is >90%.

[0254] The chessboard has a period of 150km, a grid resolution of about 50km, a spatial positioning error of <10%, and a periodic feature recovery rate of >90%.

[0255] Computational efficiency: The total computation time for 32 MPI processes is 12 hours, which is suitable for efficient inversion of large-scale models.

[0256] Verification Case 2: Inversion of Localized Blocky Anomalous Bodies

[0257] Model and data preparation:

[0258] Basic configuration: Uses the same spherical shell mesh (131,842 elements) and excitation source and station configuration (3 sources) as the chessboard model. (1200 stations per second).

[0259] Conductivity settings: Background conductivity: S / m, such as Figure 5 As shown, the anomaly is a localized blocky structure located at longitude 20°–60°W, latitude 30°–60°N, and depth 100–300 km. The anomaly exhibits a concentrated distribution relative to the background and has an electrical conductivity 10 times that of the background. S / m.

[0260] Observational data: Same as the chessboard model, with 1% Gaussian noise added.

[0261] Inversion parameter settings:

[0262] Based on the characteristics of local anomalies, we use differentiated inversion parameters compared to the first verification case to suppress false anomalies in the background and improve the accuracy of anomaly localization.

[0263] Regularization: (Strong regularization to prevent false anomalies in the background), decay rate 0.7 (slow decay).

[0264] Target RMS error: 0.5 (relatively lenient, prioritizing the recovery of anomalous features).

[0265] Computing resources: 128 MPI processes (to improve the efficiency of local detail calculations).

[0266] Inversion results and analysis:

[0267] Model feature recovery: Figure 6 The inversion evolution process shown clearly demonstrates the gradual identification and refinement of local blocky anomalies. The anomalies become clearly visible after 50-100 iterations, and are gradually refined in later stages. The main location of the anomaly is consistent with the true value, with a spatial positioning error of <500km. The conductivity of the inverted anomaly is 8-12 times that of the background, with an intensity recovery accuracy of ±20%. Small spurious anomalies in the early iterations are effectively suppressed during the regularization decay process, and the background region remains approximately uniform. Snapshots from multiple iterations demonstrate the spatial positioning capability and background suppression ability.

[0268] Iterative convergence characteristics: The iterative process (summary) is shown in Table 2. As shown in Table 2, the RMS error monotonically decreased from the initial 482.9 to 0.053, and the objective function... The objective function shows a stable decreasing trend as it decreases from 3.02 × 10¹⁰ to 2.07 × 10⁴, as shown below. Figure 7 As shown, the convergence curve reflects the inversion convergence law of local blocky anomalies. Strong initial regularization effectively avoids model divergence, and the step size is gradually adjusted with iteration, resulting in a stable inversion process. Compared with the chessboard model, the convergence speed and final accuracy are slightly different, reflecting the differences in inversion difficulty for different anomaly types.

[0269] Table 2. Iterative Process of Verification Case 2

[0270] ;

[0271] Background noise suppression: Minor spurious anomalies generated in the early iterations are effectively suppressed during the gradual decay of regularization, and the background region eventually remains approximately uniform, demonstrating good model stability.

[0272] To verify the accuracy of the adjoint gradient method, a comparative verification was performed using the finite difference method: At a specified step in the inversion iteration (e.g., the 10th iteration of the chessboard model), a subset of parameters (e.g., 50) were selected, and the change in the objective function was calculated by perturbing the parameters in small increments, thus obtaining the finite difference gradient. The gradient obtained by the adjoint gradient method In comparison, the formula is:

[0273] ;

[0274] ;

[0275] In the formula, Let be the finite difference gradient of the i-th parameter. Let m be the inversion objective function, and m be the unconstrained parameter vector. For the small numerical quantity used for parameter perturbation, This is a standard unit basis vector where only the i-th component is 1 and the rest are 0, used to individually perturb the i-th parameter; Let be the adjoint gradient of the i-th parameter. The electromagnetic field angular frequency, The permeability of free space, For extracting the real part, Let be the conjugate transpose of the adjoint field vector corresponding to the i-th parameter. Let be the partial derivative of the mass matrix with respect to the i-th parameter. This is the forward electric field coefficient vector, and the negative sign is used to ensure that the gradient follows the descent direction of the objective function.

[0276] The verification results show that the average relative error of the gradient is 0.3%, and the errors of all parameters are all <5%, confirming the accuracy of the accompanying gradient calculation and providing reliable gradient information for subsequent optimization iterations.

[0277] Example 2:

[0278] like Figure 8 As shown, this embodiment provides a spherical shell three-dimensional magnetotelluric conductivity inversion system for implementing the method of Embodiment 1, including an input module, a forward calculation module, an inversion calculation module, an optimization module, and an output module;

[0279] The input module is used to perform the following: constructing a three-dimensional finite element mesh for the spherical shell, setting the initial conductivity value of each element in the three-dimensional finite element mesh for the spherical shell, establishing the mapping relationship between the unconstrained parameter vector and the conductivity, and obtaining the initial conductivity model of the spherical shell region;

[0280] The forward modeling module is used to perform the following: based on the conductivity model, construct the corresponding frequency domain Maxwell equations; based on the frequency domain Maxwell equations, perform finite element discrete frequency domain forward modeling for each excitation source and frequency, and solve for the electric field coefficient vector.

[0281] Based on the three-dimensional finite element mesh of the spherical shell and the electric field coefficient vector, three-dimensional spatial interpolation is performed at the location of the observation station, and after coordinate transformation and Faraday's law derivation, the multi-source tangential electromagnetic field components of each station are extracted.

[0282] The inversion calculation module is used to perform the following: based on the multi-source tangential electromagnetic field components of each station, the least squares formula is used to calculate and construct the multi-source impedance tensor.

[0283] Based on the multi-source impedance tensor, a total objective function composed of a weighted data fitting term and a regularization term is defined;

[0284] Based on the derivative of the data fitting term with respect to the electric field coefficient vector, the adjoint equation is constructed and the adjoint field vector is obtained by solving it; combined with the adjoint field vector, the partial derivative of the objective function with respect to the parameter vector is calculated.

[0285] The optimization module is used to perform the following: calculate the search direction and determine the iteration step size based on the partial derivative of the objective function with respect to the parameter vector, and update the parameter vector; perform iterative calculations through the forward calculation module and the inverse calculation module until the convergence condition is met or the maximum number of iterations is reached;

[0286] The output module is used to output the final conductivity model of the spherical shell region.

[0287] In one possible implementation, the system further includes a parallel module for implementing parallel computing.

[0288] In the above scheme, the system reads basic data through the input module, performs frequency domain forward modeling and electromagnetic field interpolation through the forward modeling module, constructs the impedance tensor and calculates the parameter gradient through the inversion modeling module, implements parameter updates and regularization adjustments through the optimization module, determines convergence and outputs the conductivity model through the output module, and implements MPI parallelism, global gradient reduction and load balancing through the parallel module.

[0289] In some embodiments, the input module may include: a spherical shell finite element mesh reader (supporting Gmsh format), an initial conductivity value reader (supporting element-by-element specification), an excitation source file parser (supporting volume current density cell and wise formats), and an observation data reader (supporting multiple frequencies, multiple polarizations, and multiple stations).

[0290] The forward modeling module may include: a finite element matrix assembler ( and Matrix, Nédélec edge element), frequency domain system matrix constructor (automatically connects to different frequencies) ), FGMRES linear solver (integrated AMS preconditions), station field interpolator (3D spatial coordinate transformation);

[0291] Inversion Calculation Module: Multi-source Impedance Tensor Calculator, Adjoint Source Constructor (Automatic Impedance Jacobian Derivation), Adjoint Equation Solver, Gradient Calculator and Accumulator;

[0292] The optimization module may include: an L-BFGS optimizer (managing historical vectors), a line search module (Armijo conditions), an adaptive regularized parameter decay controller, and a parameterized to non-parameterized conversion interface;

[0293] The output module may include: conductivity model output, iterative summary recorder, and intermediate checkpoint saver (supporting breakpoint restart);

[0294] Parallel modules may include: MPI communication interface, global gradient optimization, and load balancer.

[0295] In some embodiments, the system further integrates a parameterization conversion module for implementing Logistic parameterization mapping and reciprocal calculation, supporting multi-segment parameterization mapping, with each segment employing an independent physical lower bound on conductivity. and the Upper Realm .

[0296] Example 3:

[0297] This embodiment provides an electronic device, including: a memory and a processor;

[0298] The memory is used to store computer programs;

[0299] The processor is configured to invoke the computer program to execute the method as described in Embodiment 1.

[0300] Example 4:

[0301] This embodiment provides a computer-readable storage medium storing a computer program. When the computer program is run on an electronic device, it causes the electronic device to perform the method described in Embodiment 1.

[0302] Example 5:

[0303] This embodiment provides a computer program product, including a computer program that, when run on an electronic device, causes the electronic device to perform the method described in Embodiment 1.

[0304] The specific implementation of the system, electronic device, computer-readable storage medium, and computer program product provided in this application can be referred to the specific embodiments of the above methods, and will not be repeated here.

[0305] Obviously, those skilled in the art should understand that the various units or steps of this application described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by a computing device, or fabricating them separately as individual integrated circuit modules, or fabricating multiple modules or steps into a single integrated circuit module. Thus, this application is not limited to any particular combination of hardware and software.

[0306] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method of 3D magnetotelluric conductivity inversion for a spherical shell, characterized by, include: S1. Construct a three-dimensional finite element mesh for the spherical shell, set the initial conductivity value of each element in the three-dimensional finite element mesh for the spherical shell, establish the mapping relationship between the unconstrained parameter vector and the conductivity, and obtain the initial conductivity model of the spherical shell region. S2. Based on the conductivity model, construct the corresponding frequency domain Maxwell equations; based on the frequency domain Maxwell equations, perform finite element discrete frequency domain forward modeling for each excitation source and frequency, and solve for the electric field coefficient vector. S3. Based on the three-dimensional finite element mesh of the spherical shell and the electric field coefficient vector, three-dimensional spatial interpolation is performed at the location of the observation station, and after coordinate transformation and Faraday's law derivation, the multi-source tangential electromagnetic field components of each station are extracted. S4. Based on the multi-source tangential electromagnetic field components of each station, the multi-source impedance tensor is constructed by calculating the least squares formula. S5. Based on the multi-source impedance tensor, define the overall objective function composed of the weighted data fitting term and the regularization term; S6. Based on the derivative of the data fitting term with respect to the electric field coefficient vector, construct the adjoint equation and solve for the adjoint field vector; combine the adjoint field vector to calculate the partial derivative of the objective function with respect to the parameter vector. S7. Based on the partial derivative of the objective function with respect to the parameter vector, calculate the search direction and determine the iteration step size, and update the parameter vector; S8. Repeat steps S2 to S7 for iterative calculations until the convergence condition is met or the maximum number of iterations is reached, and output the final conductivity model of the spherical shell region.

2. The method according to claim 1, characterized in that, In step S1, the mapping relationship between the unconstrained parameter vector and conductivity is established using Logistic parameterization, and the expression is: ； The derivative is: ; in, This represents the conductivity of the r-th grid cell. and These are the lower and upper physical bounds of conductivity, respectively. It is a parameter vector The first in One parameter.

3. The method according to claim 1, characterized in that, In step S2, when performing finite element discrete frequency-domain forward modeling of each excitation source and frequency based on the frequency-domain Maxwell equations, the specific steps include: discretizing the frequency-domain Maxwell equations into a complex linear system using finite element methods. ,in For the system matrix corresponding to the frequency, Let be the electric field coefficient vector corresponding to the s-th excitation source. This is the source term vector corresponding to the s-th excitation source; The complex linear system is solved using the preconditioned FGMRES iterative solver to obtain the electric field coefficient vector; For the matrix and triangular decomposition results of a multi-excitation source multiplexing system with the same frequency, only the source term vector corresponding to each excitation source is updated; MPI parallel computing is used, with each process handling different excitation sources or grid subdomains.

4. The method according to claim 1, characterized in that, In step S3, the specific process of performing three-dimensional spatial interpolation at the observation station location and extracting the multi-source tangential electromagnetic field components of each station through coordinate transformation and Faraday's law derivation is as follows: calculate the centroid coordinates of the tetrahedral element where the station is located, calculate the Cartesian component of the electromagnetic field at the station location using the Nédélec edge element basis function, convert it to the tangential electric field component through the spherical coordinate Jacobian matrix, and then derive the tangential magnetic field component from the electric field curl according to Faraday's law.

5. The method according to claim 1, characterized in that, The least squares formula for constructing the multi-source impedance tensor in step S4 is as follows: ； in, Let k represent the magnetotelluric impedance tensor of station k. Indicates conjugate transpose; and The tangential electric field matrix and tangential magnetic field matrix are composed of the multi-source tangential electric field and magnetic field components of station k.

6. The method according to claim 1, characterized in that, The overall objective function in step S5 is: ； ； ； in, Let be the overall objective function. For the first The regularization parameter for the next iteration; This represents the objective function for fitting the data. Represents the parameter vector The corresponding conductivity model, the predicted impedance real vector obtained through forward modeling and impedance construction; Represents the measured real vector of impedance; Represents a data weighting matrix; This represents the standard function for regularization. This represents the three-dimensional solution domain of the spherical shell. Indicates spatial location Depth weight at the location, Represents the parameter vector In spatial location The gradient at that point.

7. The method according to claim 1, characterized in that, Step S6 specifically includes: S61. Link the derivatives of the data fitting terms with respect to the impedance tensor to the derivatives with respect to the electric field coefficient vector using the chain rule, construct the adjoint equation, and solve it to obtain the adjoint field vectors corresponding to each excitation source. S62. Combine the partial derivatives of the adjoint field vector and the mass matrix with respect to the parameter vector, and calculate the partial derivatives of the data fitting term with respect to the parameter vector using the chain rule. S63. Calculate the gradient of the regularization term with respect to the parameter vector, and sum the gradient of the data fitting term and the gradient of the regularization term with weights to obtain the partial derivative of the total objective function with respect to the parameter vector. The adjoint equation is in the form of: ; In the formula, Represents the conjugate transpose of the system matrix. This represents the adjoint field vector corresponding to the s-th excitation source. This represents the electric field coefficient vector corresponding to the s-th excitation source. This is the partial derivative of the data fitting term with respect to the electric field coefficient vector; The data fitting term applies to the r-th parameter in the parameter vector. The partial derivatives are calculated using the adjoint-state method, and the formula is: ； In the formula, The electromagnetic field angular frequency, The permeability of free space, Operator for extracting the real part; Let mr be the partial derivative of the mass matrix with respect to the r-th parameter, and this partial derivative is non-zero only within the mesh cell corresponding to mr.

8. The method according to claim 1, characterized in that, In step S7, the search direction is calculated using the L-BFGS algorithm, and the iteration step size is determined using the Armijo line search.

9. The method according to claim 1, characterized in that, During the iteration of S8, the regularization parameter of the regularization term in step S5 is dynamically adjusted based on the decrease in RMS error after the parameter vector update. The dynamic adjustment employs an exponential decay strategy, expressed as follows: ； in, Let be the regularization parameter for the nth iteration, where n is the iteration number. These are the initial values ​​for the regularization parameters. This is the lower bound of the regularization parameter. The decay rate is the regularization parameter.

10. A three-dimensional magnetotelluric conductivity inversion system for a spherical shell, characterized in that, The method for performing any one of the methods of claims 1-9 includes an input module, a forward calculation module, an inverse calculation module, an optimization module, and an output module; The input module is used to perform the following: constructing a three-dimensional finite element mesh for the spherical shell, setting the initial conductivity value of each element in the three-dimensional finite element mesh for the spherical shell, establishing the mapping relationship between the unconstrained parameter vector and the conductivity, and obtaining the initial conductivity model of the spherical shell region; The forward modeling module is used to perform the following: constructing the corresponding frequency-domain Maxwell equations based on the conductivity model; Based on Maxwell's equations in the frequency domain, finite element discrete frequency domain forward modeling is performed on each excitation source and frequency to obtain the electric field coefficient vector. Based on the three-dimensional finite element mesh of the spherical shell and the electric field coefficient vector, three-dimensional spatial interpolation is performed at the location of the observation station, and after coordinate transformation and Faraday's law derivation, the multi-source tangential electromagnetic field components of each station are extracted. The inversion calculation module is used to perform the following: based on the multi-source tangential electromagnetic field components of each station, the least squares formula is used to calculate and construct the multi-source impedance tensor. Based on the multi-source impedance tensor, a total objective function composed of a weighted data fitting term and a regularization term is defined; Based on the derivative of the data fitting term with respect to the electric field coefficient vector, the adjoint equation is constructed and the adjoint field vector is obtained by solving it; combined with the adjoint field vector, the partial derivative of the objective function with respect to the parameter vector is calculated. The optimization module is used to perform the following: calculate the search direction and determine the iteration step size based on the partial derivative of the objective function with respect to the parameter vector, and update the parameter vector; Iterative calculations are performed using the forward and inverse calculation modules until the convergence condition is met or the maximum number of iterations is reached. The output module is used to output the final conductivity model of the spherical shell region.

Images