CN-122021210-A - Spherical shell three-dimensional magnetotelluric conductivity inversion method and system

CN122021210ACN 122021210 ACN122021210 ACN 122021210ACN-122021210-A

Abstract

The invention discloses a spherical shell three-dimensional magnetotelluric conductivity inversion method which comprises the steps of constructing a spherical shell three-dimensional finite element grid, establishing a mapping relation between unconstrained parameter vectors and conductivity to obtain an initial conductivity model, constructing a frequency domain Maxwell equation based on the model, performing frequency domain forward modeling on excitation sources and frequencies after finite element dispersion to solve electric field coefficient vectors, interpolating at an observation station position, deriving through coordinate transformation and Faraday law to extract multisource tangential electromagnetic field components, constructing a multisource impedance tensor by adopting a least square formula, defining a total objective function containing weighted data fitting items and regularization items, solving accompanying field vectors through an accompanying equation, calculating partial derivatives of the objective function on the parameter vectors, determining a search direction and iteration step length based on the partial derivatives, updating the parameter vectors, and outputting a spherical shell region final conductivity model after iteration to convergence. The method breaks through the limitation of traditional plane inversion, adapts to full sphere scale application, and improves inversion efficiency and accuracy.

Inventors

Wang Juanyu
CHEN CHAOJIAN
LI RENJIE
XIE LIANGYU
Xiang Jingnian

Assignees

中南大学

Dates

Publication Date: 20260512
Application Date: 20260416

Claims (10)

1. The spherical shell three-dimensional magnetotelluric conductivity inversion method is characterized by comprising the following steps of: S1, constructing a spherical shell three-dimensional finite element grid, setting initial conductivity values of all units in the spherical shell three-dimensional finite element grid, and establishing a mapping relation between unconstrained parameter vectors and conductivity to obtain an initial conductivity model of a spherical shell region; S2, constructing a corresponding frequency domain Maxwell equation based on the conductivity model, performing finite element discrete frequency forward modeling on each excitation source and frequency based on the frequency domain Maxwell equation, and solving to obtain an electric field coefficient vector; S3, based on spherical shell three-dimensional finite element grids and electric field coefficient vectors, performing three-dimensional spatial interpolation at the positions of observation stations, and extracting multisource tangential electromagnetic field components of each station through coordinate transformation and Faraday law derivation; s4, calculating and constructing a multisource impedance tensor by adopting a least square formula based on multisource tangential electromagnetic field components of each station; s5, defining a total objective function consisting of a weighted data fitting term and a regularization term based on the multi-source impedance tensor; s6, constructing an accompanying equation based on the derivative of the data fitting term on the electric field coefficient vector, and solving to obtain an accompanying field vector; S7, calculating a search direction and determining an iteration step length based on the partial derivative of the objective function on the parameter vector, and updating the parameter vector; and S8, repeating the steps S2 to S7 to perform iterative computation until convergence conditions are met or the maximum iteration times are reached, and outputting a final conductivity model of the spherical shell region.
2. The method of claim 1, wherein establishing the unconstrained mapping relationship between the parameter vector and the conductivity in step S1 uses Logistic parameterized mapping, expressed as: ; The derivative is: ; Wherein, the Indicating the conductivity of the nth grid cell, And The physical lower and upper boundaries of conductivity are respectively, Is a parameter vector The first of (3) And parameters.
3. The method of claim 1, wherein in step S2, finite element discrete frequency forward modeling is performed on each excitation source and frequency based on frequency domain Maxwell 'S equations, specifically comprising the steps of discretizing the frequency domain Maxwell' S equations into a complex linear system through finite elements Wherein For a system matrix of corresponding frequencies, For the electric field coefficient vector corresponding to the s-th excitation source, The source item vector corresponding to the s-th excitation source; solving the complex linear system through a precondition FGMRES iteration solver to obtain an electric field coefficient vector; multiplexing the system matrix and the triangular decomposition result by multiple excitation sources with the same frequency, and only updating the source item vectors corresponding to the excitation sources; MPI parallel computation is adopted, and each process processes different excitation sources or grid subdomains.
4. The method of claim 1, wherein in step S3, three-dimensional spatial interpolation is performed at the position of the observation station and the tangential electromagnetic field component of each station is extracted through coordinate transformation and Faraday ' S law, and the specific process of extracting the tangential electromagnetic field component of each station is that the barycentric coordinates of the tetrahedral unit where the station is located are calculated, the Cartesian component of the electromagnetic field at the position of the station is calculated by using the Lec ' S primitive function, the Cartesian component is converted into the tangential electric field component through the spherical coordinate Jacobian matrix, and the tangential magnetic field component is derived from the electric field rotation according to Faraday ' S law.
5. The method of claim 1, wherein constructing a least squares formulation for the multisource impedance tensor in step S4 is: ; Wherein, the Representing the magnetotelluric impedance tensor of station number k, Represents a conjugate transpose; And A tangential electric field matrix and a tangential magnetic field matrix which are formed by multisource tangential electric fields and magnetic field components of the k-number station.
6. The method according to claim 1, wherein the total objective function in step S5 is: ; ; ; Wherein, the As a function of the overall objective function, Is the first Regularization parameters of the secondary iterations; Representing a data fitting project label function; Representing the result of the parameter vector A corresponding conductivity model, a predicted impedance real vector obtained through forward modeling and impedance construction; Representing an actual measured impedance vector; representing a data weighting matrix; Representing the regularized project label function, Representing the sphere three-dimensional solving domain, Representing spatial position The depth weight at which to place, Representing parameter vectors In a space position A gradient thereat.
7. The method according to claim 1, wherein step S6 specifically comprises: S61, the derivative of the data fitting item on the impedance tensor is linked to the derivative of the electric field coefficient vector through a chain rule, an accompanying equation is constructed, and the accompanying field vector corresponding to each excitation source is obtained through solving; S62, combining the partial derivatives of the parameter vectors by the accompanying field vectors and the quality matrix, and calculating the partial derivatives of the parameter vectors by the data fitting items through a chain rule; S63, calculating the gradient of the regularization term on the parameter vector, and carrying out weighted summation on the gradient of the data fitting term and the regularization term to obtain the partial derivative of the total objective function on the parameter vector; Wherein the accompanying equation is in the form of ; In the formula, Representing the conjugate transpose of the system matrix, Representing the companion field vector corresponding to the s-th excitation source, Represents the electric field coefficient vector corresponding to the s-th excitation source, Partial derivatives of the data fitting term to the electric field coefficient vector; the data fitting term pairs the r-th parameter in the parameter vector The partial derivative of (2) is calculated by the adjoint-state method, and the formula is: ; In the formula, For the angular frequency of the electromagnetic field, Is the magnetic permeability of the vacuum and is equal to the magnetic permeability of the vacuum, Extracting an operator for a real part; the partial derivative of the quality matrix with respect to the r-th parameter is non-zero only within the grid cell to which mr corresponds.
8. The method according to claim 1, wherein in step S7, the search direction is calculated by the L-BFGS algorithm, and the iteration step is determined by Armijo line search.
9. The method of claim 1, wherein the regularization parameters of the regularization term in step S5 are dynamically adjusted based on the magnitude of the RMS error decrease after the parameter vector update during the iteration of S8 Wherein the dynamic adjustment adopts exponential decay strategy adjustment, and the expression is: ; Wherein, the Regularization parameter of the nth iteration, n is the iteration number, For the initial value of the regularization parameter, For the regularized lower bound of the parameters, Is the regularized parameter decay rate.
10. A spherical shell three-dimensional magnetotelluric conductivity inversion system, which is characterized by being used for executing the method of any one of claims 1-9, and comprising an input module, a forward calculation module, an inversion calculation module, an optimization module and an output module; the input module is used for executing the steps of constructing a spherical shell three-dimensional finite element grid, setting initial conductivity values of all units in the spherical shell three-dimensional finite element grid, and establishing a mapping relation between unconstrained parameter vectors and conductivity to obtain an initial conductivity model of a spherical shell region; The forward calculation module is used for executing the finite element discrete frequency domain forward calculation of each excitation source and frequency based on the frequency domain Maxwell equation to obtain an electric field coefficient vector by solving; Based on spherical shell three-dimensional finite element grids and electric field coefficient vectors, three-dimensional spatial interpolation is carried out at the positions of observation stations, and multi-source tangential electromagnetic field components of each station are extracted through coordinate transformation and Faraday law derivation; the inversion calculation module is used for executing the calculation and construction of a multi-source impedance tensor by adopting a least square formula based on multi-source tangential electromagnetic field components of each station; defining a total objective function consisting of a weighted data fitting term and a regularization term based on the multi-source impedance tensor; constructing an accompanying equation based on the derivative of the data fitting term on the electric field coefficient vector, and solving to obtain an accompanying field vector; the optimization module is used for executing the iterative computation through the forward computation module and the inversion computation module until the convergence condition is met or the maximum iteration times are reached; And the output module is used for outputting the final conductivity model of the spherical shell region.

Description

Spherical shell three-dimensional magnetotelluric conductivity inversion method and system Technical Field The invention relates to the technical field of geophysical inversion, in particular to a three-dimensional inversion method of a Magnetotelluric (MT). Specifically, the conductivity distribution of the spherical shell region is inverted by finite element forward and concomitant gradient methods using multisource body current excitation and multi-station impedance observation data. Background Magnetotelluric (Magnetotelluric, MT) is a passive electromagnetic detection method. Researchers reversely push out the conductivity structure of the underground medium by recording the change of the surface natural electromagnetic field along with time. In other words, we do not actively transmit signals, but rather use electromagnetic disturbances of the earth itself to obtain subsurface information. In practical applications, the range of use of MT is very wide. From basic science to resource exploration, it plays an important role. For example, in global scale studies, researchers use MT data to infer mantle temperature distribution, water content, and the presence of a portion of the melt, MT is commonly used in the field of engineering geophysics for the detection of mineral resources and geothermal systems, in structural geology studies, it can reveal plate boundaries, fracture zones, and deep structural features, and in hydrocarbon exploration, conductivity structures are also commonly used to identify electrical anomalies in deep reservoirs. MT inversion software, such as OCCAM, rhoplus and ModEM, which are widely used today, are not quite different in basic idea. Most systems assume that the incident electromagnetic field is a plane wave and establish the electromagnetic field equation under a cartesian coordinate system. The computational object is typically a horizontal hierarchical model, or a two-dimensional, three-dimensional conductivity structure at the region scale. To simplify the problem, these methods often employ single source, single polarization excitation assumptions, and solve model parameters by Tikhonov regularization in combination with damped least squares iteration. Within the area scale (typically less than about 500 km), such methods have been validated by a number of practices and have formed a relatively mature technological system. The above approach gradually exposes significant limitations as the scope of research extends to the global scale. First is the problem of coordinate system and geometric description. Conventional MT inversion is typically based on plane wave assumptions and modeled in a cartesian coordinate system. But at the global scale the effect of the earth's curvature cannot be neglected anymore. The observation stations are distributed on the spherical surface, and the underground abnormal body also exists in the spherical shell structure. If planar approximations are still used, it is difficult to accurately describe the true geometric relationships, and such methods tend to produce significant errors in processing data across a large geographic range. Another problem comes from the electromagnetic excitation itself. The global electromagnetic field is not a simple single-source single-polarization structure, but is a complex system formed by superposition of multiple sources, such as magnetic layer current, ion layer current, tidal magnetic field, and the like. In this case, the applicability of the conventional formula (z=e/H, where Z is the scalar impedance, E is the horizontal tangential electric field scalar value observed at the earth's surface, and H is the horizontal tangential magnetic field scalar value orthogonal to the electric field) is limited, since the formula implies the assumption of a single excitation source. For multi-source, multi-polarized electromagnetic fields, there is still a lack of uniform and clear processing methods how the impedance tensor should be defined and how to stabilize the calculation. This makes it difficult for the prior art to fully utilize all of the information contained in the observed data during the actual inversion process. The problem is also apparent in terms of parameter sensitivity calculations. If inversion is performed in the spherical shell coordinate system, the observation station does not always correspond to the node position in the finite element or finite difference grid, and thus it is necessary to obtain electric field and magnetic field values by three-dimensional spatial interpolation. Meanwhile, the derivative of the multi-source impedance tensor on the conductivity parameter relates to a complex chained derivation process, and the calculation steps are complex. If finite difference methods are employed to estimate sensitivity one by one perturbing model parameters, a large number of forward calculations often need to be performed. In a global scale model, this computational effort is often u