A magnetotelluric simulation method for an arbitrary anisotropic conductivity tensor in spherical coordinates

Abstract

The present invention discloses a magnetotelluric simulation method for an arbitrary anisotropic conductivity tensor in spherical coordinate system, belonging to the technical field of geophysical forward modeling, including: determining the coordinates and observation frequencies of a measuring point array; constructing an earth background model and discretizing it using unstructured tetrahedral meshes to obtain unstructured meshes; applying an arbitrary anisotropic conductivity tensor in spherical coordinates to assign electrical properties to the unstructured meshes; applying orthogonal non-planar wave field sources to the outer meshes of the unstructured meshes and calculating the source term integrals of the non-planar wave field sources; calculating the mass matrix and stiffness matrix of the elements in the unstructured meshes and forming the overall element matrix according to the numbering relationship; constructing a corresponding simulation equation set according to the overall element matrix and the source term integrals and at the current observation frequency; solving the simulation equation set to achieve magnetotelluric simulation, so as to solve the technical problem that the existing magnetotelluric simulation technology cannot flexibly simulate the magnetotelluric responses of arbitrarily anisotropic media at the intercontinental scale.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical forward modeling, and particularly relates to a magnetotelluric simulation method for arbitrary anisotropic conductivity tensors in spherical coordinates. Background Technique

[0002] The activities and evolution of the Earth not only bring resources to humans but also accompany natural disasters. Studying scientific issues related to the Earth is of great significance for the sustainable development of humanity. Geophysics is an effective method for detecting the distribution of the Earth's material structure through physical fields and has been widely used in disaster monitoring, resource exploration, and the study of Earth science issues. The detection methods in the geophysical field for the mantle scale mainly include seismic observation and magnetotelluric sounding. Among them, the magnetotelluric sounding method uses natural field sources for underground detection and has many advantages such as large exploration depth and low cost, and is a sharp tool for deep Earth exploration. Nowadays, in order to fully utilize the detection advantages of the magnetotelluric method and achieve accurate imaging of the electrical structure in the deep Earth, a number of intercontinental-scale magnetotelluric sounding plans have been deployed internationally. Intercontinental-scale magnetotelluric observations can achieve the overall imaging of underground media on a large scale and reduce non-uniqueness, but at the same time, they also bring new challenges to data inversion.

[0003] Forward simulation is the basis of inversion. As the scale of magnetotelluric observations becomes larger and larger, the magnetotelluric simulation in spherical coordinates has become a hot and difficult point in the field of geophysical research. First of all, large-scale simulations often require discretized solutions, and the spherical coordinate forward modeling based on the finite difference method of staggered grids has developed rapidly. Many scholars have analyzed the differences in the magnetotelluric responses in spherical coordinates and Cartesian coordinates and have unanimously agreed that the influence of the Earth's curvature on large-scale long-period magnetotelluric data cannot be ignored. The finite element method based on octree grids has realized the magnetotelluric forward modeling in spherical coordinates. However, considering that complex land-sea boundaries, the Earth's core, and the ionosphere have a significant impact on the electromagnetic response, global-scale magnetotelluric forward modeling requires multi-scale model discretization techniques.

[0004] The unstructured tetrahedron can flexibly fit structures of various scales and is an optimal method for global magnetotelluric 3D forward modeling. Secondly, as the simulation scale increases, the form of the field source changes significantly. The traditional magnetotelluric field source is assumed to come from the ionospheric system far from the Earth's surface (at 4 - 10 times the Earth's radius). For small-scale survey areas, it can be considered that a uniform plane polarized wave is vertically incident on the Earth's medium. However, when the survey area is of intercontinental scale, this assumption will no longer be applicable, and a spherical coordinate system field source conforming to the ionospheric structure should be established at this time. Based on the finite difference algorithm, magnetotelluric simulation of isotropic media in the spherical coordinate system is realized, and the response differences between the plane wave field source and the non-uniform electric field source are compared, believing that the field source will affect long-period data. Finally, the inhomogeneity of the mantle leads to electrical anisotropy, which will significantly affect the magnetotelluric observation signal. The existing magnetotelluric methods in the spherical coordinate system are based on the isotropic assumption and are prone to bring incorrect results. There is less research on simulating the magnetotelluric response of electrically anisotropic media in the spherical coordinate system. One is to use the conductivity anisotropy tensor in the Cartesian coordinate system in the spherical coordinate system to simulate the geomagnetic sounding response, but this obviously has serious deviations, and the elements of its anisotropy tensor are discontinuous along the longitude and latitude directions.

[0005] Therefore, in order to accurately simulate the magnetotelluric field in the spherical coordinate system to meet the needs of processing intercontinental magnetotelluric array data in deep Earth observations, it is urgent to develop a 3D simulation method for magnetotellurics in arbitrarily anisotropic media in the spherical coordinate system to serve the inversion of large-scale array data and deep Earth scientific research in the future. Summary of the Invention

[0006] Aiming at the technical problem that the existing magnetotelluric simulation technology cannot flexibly simulate the magnetotelluric response of arbitrarily anisotropic media at the intercontinental scale, the present invention proposes a magnetotelluric simulation method for an arbitrarily anisotropic conductivity tensor in the spherical coordinate system. It uses unstructured tetrahedral meshes to flexibly depict the real Earth environment and uses spherical harmonic functions to simulate the magnetotelluric orthogonal field source, proposes a new conductivity tensor for arbitrarily anisotropic media in the spherical coordinate system and its equivalent form in the Cartesian coordinate system, and finally gives the simulation equation and the solution method based on the vector finite element method. It provides a forward modeling technology for the inversion research of intercontinental scale array magnetotelluric data and is expected to serve deep Earth scientific research in the future.

[0007] To achieve the above-mentioned invention purpose, an embodiment provides a magnetotelluric simulation method for an arbitrarily anisotropic conductivity tensor in the spherical coordinate system, including the following steps:

[0008] Determine the coordinates and observation frequencies of the measurement point array according to the magnetotelluric simulation target;

[0009] Construct an Earth background model and discretize it using unstructured tetrahedral meshes to obtain unstructured meshes;

[0010] Assign electrical properties to unstructured grids using an arbitrary anisotropic conductivity tensor in spherical coordinates;

[0011] Apply an orthogonal non-planar wave field source to the outer grids of the unstructured grid and calculate the source term integral of the non-planar wave field source;

[0012] Calculate the mass matrix and stiffness matrix of the elements in the unstructured grid and form the overall element matrix according to the numbering relationship;

[0013] Construct a corresponding simulation equation set according to the overall element matrix and the source term integral and at the current observation frequency;

[0014] Solve the simulation equation set to obtain the electric fields of all element edges and calculate the magnetotelluric response data based on the electric fields through interpolation functions to achieve magnetotelluric simulation.

[0015] Preferably, assigning electrical properties to unstructured grids using an arbitrary anisotropic conductivity tensor in spherical coordinates includes:

[0016] Give its expression in spherical coordinates as:

[0017] ;

[0018] Where the superscript s Indicates that the physical quantity is in spherical coordinates, that is Indicates the anisotropic conductivity tensor in spherical coordinates, and the right side of the formula is The expansion form of, where the subscript of each component indicates providing its corresponding scalar conductivity along different spherical coordinate system directions, and the tensor will be further converted into three spherical coordinate principal axis conductivities using spatial Euler rotation And three spherical coordinate principal axis rotation angles To describe:

[0019] ;

[0020] Where the superscript T represents the transpose of the matrix, and , , Are the three principal axes in spherical coordinates , And The corresponding spatial Euler rotation matrices, expressed as:

[0021] , , ;

[0022] The arbitrary anisotropic conductivity tensor in the spherical coordinate system satisfies Ohm's law , where the superscript s Indicates that the physical quantity is in spherical coordinates, represents the current density in spherical coordinates, represents the electric field in spherical coordinates;

[0023] gives the Cartesian coordinate equivalent form of the arbitrary anisotropic conductivity tensor in the spherical coordinate system , where the superscript c represents the physical quantity in the Cartesian coordinate system, represents the current density in Cartesian coordinates, represents the electric field in Cartesian coordinates, and the Cartesian equivalent form of the arbitrary anisotropic conductivity tensor in the spherical coordinate system is obtained as:

[0024] ;

[0025] where represents the anisotropic conductivity tensor in the Cartesian coordinate system, and are two-direction rotation matrices, and their expressions are as follows:

[0026] , ;

[0027] Based on this, the electrical properties of the assigned unstructured grid include .

[0028] Preferably, an orthogonal non-planar wave field source is applied to the outer grid of the unstructured grid, including:

[0029] Two types of orthogonal magnetic fields are given H as the field source of the external boundary. Assuming air insulation, the magnetic field is represented by the negative gradient of the magnetic scalar potential U as:

[0030] ;

[0031] The superscript ext represents the physical quantity on the outer boundary of the grid, that is represents the magnetic field on the outer boundary, represents the magnetic scalar potential on the outer boundary, represents taking the gradient of the physical quantity, is the magnetic permeability in vacuum, and the magnetic scalar potential U in the spherical coordinate system is expressed as spherical harmonics:

[0032] ;

[0033] where represents the imaginary unit, is the radius of the earth, represents the external source coefficient, and exp represents the exponential function calculation, PFor the associated Legendre polynomials, where m and n are the degrees and orders of the polynomials respectively, in the present invention, the field distribution in the quiet solar state is taken, i.e., m = 0 and n = 1. When r is taken as 5 times the Earth's radius, the first type of field source at the outer boundary is:

[0034] ;

[0035] where, and and are the magnetic fields in three directions in the spherical coordinate system. According to the conversion relationship between spherical coordinates and Cartesian coordinates, the first type of field source in the Cartesian coordinate system ( , , ) is:

[0036] ;

[0037] The second type of field source ( , , ) takes the orthogonal field of the first type of field source as:

[0038] ;

[0039] where, and and are the magnetic fields in three directions of the first type of field source in the Cartesian coordinate system, and and are the magnetic fields in three directions of the second type of field source in the Cartesian coordinate system.

[0040] Preferably, calculating the source term integral of the non-planar wave field source includes:

[0041] Using Gaussian integration to calculate the integral of the right-hand source term. The integral term is regarded as a function related only to the spatial position and is expressed as F. Applying the integration nodes of the spatial triangle, the source term integral is:

[0042] ;

[0043] where, the superscripts 1 and 2 represent the types of field sources, i.e., and represent the source term integrals corresponding to the first type of field source and the second type of field source respectively, A is the weighting coefficient, and the subscript k represents the integration node index.

[0044] Preferably, calculating the mass matrix and stiffness matrix of the elements in the unstructured grid and forming the overall element matrix according to the numbering relationship includes:

[0045] The electric field double curl equation satisfied by the magnetotelluric field is analyzed using the vector finite element method to obtain the mass matrix of element ele and the stiffness matrix as follows:

[0046] ;

[0047] ;

[0048] where the subscript i,j represents the local number of the edge in element ele, v represents the volume of element ele, represents the j th edge corresponding to the i th edge of the vector interpolation basis function, represents the vector interpolation basis function of the j th edge, × represents taking the curl of the vector interpolation basis function, represents the arbitrary anisotropic conductivity tensor in spherical coordinates. Each edge of the element has a global number, thus forming the overall element matrices S and M.

[0049] Preferably, according to the overall element matrix and the source term integral, a corresponding simulation equation system is constructed according to the current observation frequency, including:

[0050] According to the current observation frequency f calculate the angular frequency , is the pi, and the corresponding complex frequency coefficient is , where represents the imaginary unit. Let the matrix , then the simulation equation system at the current frequency is calculated as follows:

[0051] ;

[0052] where S is the overall stiffness matrix in the overall element matrix, M is the overall mass matrix, the elements in the vector b in the right - hand side are the source term integral, and only the positions corresponding to the outer boundary edges in the vector b have values while the rest are 0. The superscripts 1, 2 represent the types of field sources, and E represents the electric field.

[0053] Preferably, solving the simulation equation system to obtain the electric fields of all element edges, including:

[0054] When the number of simulation elements is small and the computing resources are sufficient, a direct solver is used for solving. Its principle is to perform LU decomposition on K of the simulation equation system, and back - substitute for different right - hand sides of the equations to obtain the electric fields E 1 、E 2 .

[0055] Preferably, solving the simulation equations to obtain the electric fields of all unit edges, including:

[0056] When the number of simulation units is large and the computing resources are limited, the iterative method is used for parallel solution. The flexible generalized minimum residual iterative method based on the auxiliary space preconditioner and the block diagonal right preconditioner is adopted. First, the complex equations presented by the original simulation equations are transformed into equivalent real forms, the superscripts of different field sources are removed, and the general form is given:

[0057] ;

[0058] where the superscript re represents the real part of the physical quantity, and the superscript im represents the imaginary part of the physical quantity. Denote the equation in the general form as K r ·E r =b r , and the block diagonal matrix P is:

[0059] ;

[0060] Multiply the inverse matrix of P on the right side of the coefficient matrix of the general form equation, then the original equation is K r P -1 ·u = b r , the vector u is an intermediate vector, and the solution process is divided into inner and outer layer loops. The inner layer uses the auxiliary space preconditioner and the conjugate gradient method to solve the vector , q is the iterative update step vector of the intermediate vector u, and the outer layer loop uses the block diagonal matrix P corresponding to the block diagonal right preconditioner and the generalized minimum residual method to solve the electric field E r .

[0061] Preferably, based on the electric field, calculating the magnetotelluric response data through the interpolation function to realize magnetotelluric simulation, including:

[0062] Calculating the magnetotelluric response data based on the electric field through the interpolation function includes the magnetotelluric tensor impedance response, apparent resistivity, and phase at the current observation frequency.

[0063] Compared with the prior art, the beneficial effects of the present invention at least include:

[0064] The present invention uses unstructured tetrahedral meshes to flexibly depict the real Earth environment and finely depict the coastline. At the same time, its characteristic of having hypotenuses enables the simulation of irregular anomalies (such as subducting plates, fold orogenic belts, mantle plumes, etc.). The present invention uses spherical harmonic functions to construct a non-uniform orthogonal field source to break the problem that traditional plane wave field sources cannot adapt to the Earth's curvature. At the same time, the present invention proposes a new conductivity tensor for arbitrary anisotropic media in spherical coordinates. The invented tensor is continuously distributed along the spherical coordinate system, solving the problem of discontinuous conductivity tensors in the traditional Cartesian coordinate system. At the same time, the invention gives the equivalent form of the tensor in Cartesian coordinates based on Euler rotation, enabling the application of the unstructured vector finite element method to construct the magnetotelluric simulation equations for arbitrary anisotropic media in spherical coordinates. Finally, the invention gives the solution method for the electromagnetic simulation equations. It can be seen from a specific embodiment provided by the present invention that the simulation technology proposed by the invention can accurately and effectively simulate the magnetotelluric response of arbitrary anisotropic media in spherical coordinates, which provides a forward modeling technology for future magnetotelluric data inversion at the intercontinental scale and deep Earth research. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0066] Figure 1 is a flowchart of the magnetotelluric simulation method for arbitrary anisotropic conductivity tensors in spherical coordinates provided by the embodiment;

[0067] Figure 2 is a schematic diagram of the Euler rotation method of the arbitrary anisotropic conductivity tensor in spherical coordinates provided by the embodiment;

[0068] Figure 3 is a schematic diagram of the transformation from spherical coordinates to Cartesian coordinates provided by the embodiment;

[0069] Figure 4 is a schematic diagram of the model profile and the unstructured tetrahedral mesh (the air layer is removed) provided by the embodiment;

[0070] Figure 5 is a comparison diagram of the magnetotelluric responses between spherical coordinates and the analytical solution in the embodiment. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0071] To make the objectives, technical solutions, and advantages of the present invention clearer, the following further details the present invention with reference to the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and do not limit the protection scope of the present invention.

[0072] As Figure 1 shown, a magnetotelluric simulation method for an arbitrary anisotropic conductivity tensor in spherical coordinates provided by the embodiment includes the following steps:

[0073] S1. Determine the coordinates of the measuring point array and the observation frequency according to the magnetotelluric simulation target.

[0074] When actually conducting magnetotelluric observations on an intercontinental scale, a measuring network composed of measuring stations will be established in the measuring area according to the magnetotelluric simulation target. Each measuring station has corresponding longitude and latitude coordinates. The long-period magnetotelluric observation period is 10S - 100000S. According to the skin depth of the observation target range (the depth when the amplitude value decays to 1 / e, where e is the natural constant approximately equal to 2.718281828459045), n frequency points (0.1Hz - 0.00001Hz) are selected within the period. Therefore, the coordinates of the measuring point array corresponding to the simulation target are the longitude and latitude coordinates corresponding to the measuring stations, and the observation frequency is the n frequency points.

[0075] S2. Construct an earth background model and discretize it using unstructured tetrahedral meshes to obtain unstructured meshes.

[0076] The real shape of the earth is a complex irregular ellipsoid. Different ellipsoid parameters can be used to construct the earth background model according to the simulation accuracy. For the sake of completeness of description, a sphere with a radius of 6371 km is used here to construct an ideal earth. At the same time, the land and sea are depicted according to the longitude and latitude of the coastline. The earth's interior is described by a layered spherical shell for the crust, mantle, core, etc. If specific earth structures need to be simulated, such as subducting plates, mantle plumes, etc., they can also be simply added to the background mesh. The exterior of the earth is considered to be an air layer, and generally a spherical shell with 5 times the earth's radius is used to depict the outer layer of the earth model. Since the technical bottom layer of mainstream 3D modeling software is Cartesian coordinates, the conversion relationship between the spherical coordinate model and the Cartesian coordinate model is established as shown in Equations 1 - 6:

[0077] (1);

[0078] (2);

[0079] (3);

[0080] (4);

[0081] (5);

[0082] (6);

[0083] where x ,y , z corresponds to Cartesian coordinates, while r is the radial distance from the spherical surface to the center of the sphere; is the co-latitude, with a range of 0° to 180°; is the longitude, with a range of 0° to 360°.

[0084] After having the Cartesian coordinate model, an unstructured tetrahedral mesh is used for discretization. The unstructured tetrahedron can flexibly encrypt the coastline, measurement points, and areas with drastic electrical property changes to improve the accuracy of the simulation response.

[0085] S3. Apply the spherical coordinate arbitrary anisotropic conductivity tensor to assign electrical properties to the unstructured mesh.

[0086] The anisotropic conductivity tensor can describe the capabilities of both anisotropic and isotropic electrical properties simultaneously. The present invention gives its expression in spherical coordinates as:

[0087] (7);

[0088] where the superscript s represents the physical quantity in spherical coordinates, that is, represents the anisotropic conductivity tensor in spherical coordinates. The right end of formula (7) is the expansion form of . The subscript of each component indicates providing its corresponding scalar conductivity along different spherical coordinate system directions. It can be seen from formula (7) that the conductivity tensor is continuous along the spherical surface, and at the same time, this tensor is positive definite and symmetric, satisfying the energy loss law. For the convenience of simulation, the present invention will further convert the tensor into three spherical coordinate principal axis conductivities (unit Siemens per meter, S / m) and three spherical coordinate principal axis rotation angles (unit °) to describe, that is:

[0089] (8);

[0090] In formula (8), the superscript T represents the transpose of the matrix, while , , are the three principal axes , and corresponding spatial Euler rotation matrices in spherical coordinates, expressed as:

[0091] , , (9);

[0092] Figure 2The Euler rotation method of the arbitrary anisotropic conductivity tensor in spherical coordinates of the present invention is described. By using internal rotation, the principal axis conductivities are respectively rotated around the axis - after the associated rotation axis - after the associated rotation axis to the new coordinate system . The corresponding electric field is rotated in the same way. After that, the current density in the new coordinate system needs to be rotated again to obtain the equivalent current density in the original spherical coordinate system . At this time, the arbitrary anisotropic conductivity tensor in spherical coordinates satisfies Ohm's law . Here, the superscript s also represents the physical quantity in spherical coordinates. represents the current density in spherical coordinates. represents the electric field in spherical coordinates.

[0093] Since the unstructured grid is established in the Cartesian coordinate system, the present invention gives the equivalent form of the arbitrary anisotropic conductivity tensor in spherical coordinates in the Cartesian coordinate system (satisfying ). Here, the superscript c represents the physical quantity in the Cartesian coordinate system. represents the current density in Cartesian coordinates. represents the electric field in Cartesian coordinates. As Figure 3 describes the transformation of the tensor from spherical coordinates to Cartesian coordinates, that is, the tensor in spherical coordinates is first rotated around the axis by angle so that the rotated axis is parallel to the z - axis, and then rotated around the new axis by angle so that the axis is parallel to the y - axis. At this time, from the orthogonal relationship, it can be known that the axis is parallel to the x - axis. The corresponding electric field is rotated in the same way. Then, the current density parallel to the Cartesian coordinate system after rotation needs to be rotated again to obtain the current density in the original spherical coordinate system. At this time, the equivalent form of the arbitrary anisotropic conductivity tensor in spherical coordinates in the Cartesian coordinate system is:

[0094] (10);

[0095] In formula (10), represents the anisotropic conductivity tensor in the Cartesian coordinate system. and are the rotation matrices in two directions, and their expressions are as follows:

[0096] , (11);

[0097] The angle in the transformation is fixed and unique for each non-structured tetrahedral element, so there is no need to specify it during modeling. At the same time, for isotropic media, its equivalent form is consistent with the spherical coordinate form. Therefore, this transformation only needs to be applied to the elements of anisotropic media during finite element calculation. The electrical property of the non-structured tetrahedral element that needs to be assigned is 。

[0098] S4. Apply an orthogonal non-planar wave field source to the outer grid of the unstructured grid and calculate the source term integral of the non-planar wave field source.

[0099] The total field on the outermost boundary of the external air layer of the unstructured grid is formed by the superposition of the primary field and the secondary field. However, since the outer layer of the model constructed in the present invention takes a spherical shell much larger than the radius of the earth, the secondary field excited inside the earth can be considered to have decayed to zero at the outer boundary, and the total field is replaced by the primary field. The source of the geomagnetic field is unknown, and its response function is usually obtained by assuming a source to describe. For the observation of anisotropic media, the tensor impedance response function needs to be calculated. Therefore, two types of orthogonal sources need to be assumed. The present invention intends to specify two types of orthogonal magnetic fields H As the field source of the external boundary, assuming air insulation, the magnetic field can be represented by the negative gradient of the magnetic scalar potential U as:

[0100] (12);

[0101] In formula (12), the superscript ext represents the physical quantity on the outer boundary of the grid, that is represents the magnetic field on the outer boundary, represents the magnetic scalar potential on the outer boundary, represents taking the gradient of the physical quantity, is the magnetic permeability in vacuum. The magnetic scalar potential U in the spherical coordinate system can be expressed as spherical harmonic functions, that is:

[0102] (13);

[0103] Among them, represents the imaginary unit, is the radius of the earth, represents the external source coefficient, exp represents the exponential function calculation, P is the associated Legendre polynomial, m and n are the degree and order of the polynomial respectively. In the present invention, the field distribution in the quiet day state is taken, that is, m = 0, n = 1. Taking r as 5 times the radius of the earth, the first type of field source on the outer boundary is:

[0104] (14);

[0105] Among them, and and are the magnetic fields in three directions in the spherical coordinate system. According to the conversion relationship between spherical coordinates and Cartesian coordinates, it can be known that:

[0106] (15);

[0107] (16);

[0108] (17);

[0109] where and and are the magnetic fields in three directions in the Cartesian coordinate system. Therefore, the first type of field source in the Cartesian coordinate system ( , , ) is:

[0110] (18);

[0111] Therefore, the second type of field source ( , , ) takes the orthogonal field of the first type of field source as:

[0112] (19);

[0113] When constructing the simulation equation system, the right-end source term is the integral of the field source on the boundary surface , where is the unit vector interpolation basis function, is the normal vector of the surface, indicates that the integration domain is only limited to the surface of the boundary. The present invention uses Gaussian integration to calculate the integral of the right-end source term. The integral term can be regarded as a function related only to the spatial position inside and is expressed as F. Applying the integration nodes of the spatial triangle, the source term integral can be obtained as:

[0114] (20);

[0115] where the superscripts 1 and 2 represent the first and second types of field sources, that is, and respectively represent the source term integrals corresponding to the first type of field source and the second type of field source. A is the weighting coefficient, and the subscript k represents the integration node index, taking values of 1, 2, 3, 4, indicating that each triangle has four Gaussian integration nodes.

[0116] S5. Calculate the mass matrix and stiffness matrix of the elements in the unstructured grid, and form the overall element matrix according to the numbering relationship.

[0117] By analyzing the electric field double curl equation satisfied by the magnetotelluric field using the vector finite element method, the mass matrix of element ele can be obtained and the stiffness matrix are as follows

[0118] (21);

[0119] (22);

[0120] where the subscript i,j represents the local number of the edge in element ele, v represents the volume of element ele, represents the j th edge corresponding to the i th edge of the vector interpolation basis function. Each tetrahedron has 6 edges, and the number of opposite edges for each edge is also 6, including itself, represents the vector interpolation basis function of the j th edge. The symbol × represents taking the curl of the vector interpolation basis function. The edges of each element have global numbers, so the overall element matrices S and M can be formed

[0121] S6. According to the overall element matrix and the source term integration, and constructing the corresponding simulation equations according to the current observation frequency

[0122] According to the current observation frequency f calculate the angular frequency , where π is the pi, and the corresponding complex frequency coefficient is , where represents the imaginary unit. Let the matrix , then the large-scale simulation equations at the current frequency are calculated as follows

[0123] (23);

[0124] where S is the overall stiffness matrix in the overall element matrix, M is the overall mass matrix, and the elements in the vector b on the right side are the source term integration. It can be seen from S4 that only the positions corresponding to the outer boundary edges in the vector b have values and the rest are 0. The superscripts 1 and 2 represent the types of field sources, and E represents the electric field

[0125] S7. Solve the simulation equations to obtain the electric fields of all element edges

[0126] To solve the large-scale simulation equations, i.e., Equation (23), the present invention proposes to flexibly adopt different solution methods according to the scale of the equations and the computing hardware conditions. When the number of simulation units is small and the computing resources are sufficient, a direct solver is used for solution. The principle is to perform LU decomposition on K of the simulation equations, and back substitution for different right-hand sides of the equations can obtain the electric fields E of all edges. 1 and E 2 .

[0127] When the number of simulation units is large and the computing resources are limited, an iterative method is used for parallel solution. The present invention adopts a flexible generalized minimum residual iterative method based on the subspace preconditioner and the block diagonal right preconditioner. To apply this method, the complex equations presented in the original formula (23) need to be transformed into an equivalent real form. Here, the superscripts of different field sources are removed, and the general form is given:

[0128] (24);

[0129] where the superscript re represents the real part of the physical quantity, and the superscript im represents the imaginary part of the physical quantity. Denote Equation (24) as K r ·E r =b r , and the block diagonal matrix P is:

[0130] (25);

[0131] Multiply the inverse matrix of P on the right side of the coefficient matrix of Equation (24), then the original equation is K r P -1 ·u = b r . The vector u is an intermediate vector. The solution process is divided into inner and outer loops. The inner loop uses the subspace preconditioner and the conjugate gradient method to solve the vector , q is the iterative update step vector of the intermediate vector u, and the outer loop uses the block diagonal matrix P corresponding to the block diagonal right preconditioner and the flexible generalized minimum residual method to solve the electric field E r .

[0132] S8. Calculate the magnetotelluric response data based on the electric field through the interpolation function to realize magnetotelluric simulation.

[0133] Solving the large-scale simulation equations obtains the electric fields E on the edges of all unstructured tetrahedral elements. The interpolation basis functions of the unit where the measurement point is located can be used to calculate the electric fields Ex, Ey, Ez of the measurement point, and the corresponding magnetic fields can be calculated through Faraday's law of electromagnetic induction ;

[0134] The electric field E and magnetic field HIt can be calculated through the electric and magnetic fields in the Cartesian coordinate system:

[0135] (26);

[0136] (27);

[0137] (28);

[0138] (29);

[0139] (30);

[0140] (31);

[0141] During magnetotelluric observations, orthogonal stations are arranged in the northeast direction. Therefore, the observation system in which the response is located is the North (abbreviated as N)-East (abbreviated as E)-Down (abbreviated as D) coordinate system, and its relationship with the spherical coordinates is , so the magnetotelluric tensor impedance response Z at the current observation frequency is calculated as:

[0142] (32);

[0143] The apparent resistivity is:

[0144] (33);

[0145] The phase is:

[0146] (34).

[0147] S9. Repeat steps S6 - S8 until all observation frequencies are calculated, and output the simulation results.

[0148] Judge whether all frequency points have been calculated. If so, output all responses as the simulation results. Otherwise, the frequency point index nf = nf + 1, and repeat steps S6 - S8 until nf = n, that is, all observation frequencies are calculated.

[0149] To illustrate the technical effects of the above magnetotelluric simulation method for an arbitrary anisotropic conductivity tensor in the spherical coordinate system, a specific experimental example is given according to Figure 1The implementation of the magnetotelluric simulation process for an arbitrarily anisotropic conductivity tensor in the spherical coordinate system as shown. First, it is confirmed that the simulation objective of the embodiment is to verify the effectiveness of the magnetotelluric response of a three-dimensional arbitrarily anisotropic medium in the spherical coordinate system and to illustrate the necessity of using the spherical coordinate system for long-period magnetotellurics at the intercontinental scale by comparing with the analytical solution of the magnetotelluric response of a one-dimensional layered anisotropic medium. Therefore, the intersection of the equator and the prime meridian is selected as the measuring station with coordinates (0°, 0°), and 10 logarithmically equally spaced frequency points are selected as the observation frequencies from 0.1 Hz to 0.00001 Hz. Then, as Figure 4 shown, a background model of the Earth is constructed (selecting an ideal sphere with the radius of the Earth being 6371 km). To compare with the one-dimensional layered response, the interior of the Earth is simply divided into several layers, including a crust 80 km thick, a mantle anisotropic layer 220 km thick, and the remaining internal structure is equivalent to a core structure with a radius of 6071 km. After conversion to Cartesian coordinates, it is discretized using unstructured tetrahedral meshes, and local refinement is performed at the measuring point positions to obtain a mesh with a total of 4,685,766 unstructured tetrahedral elements. Then, electrical properties are assigned to the mesh. Assuming that the mantle has an anisotropic layer 220 km thick due to inhomogeneity, the corresponding spherical coordinate anisotropic parameters for this layer are set according to the tensor of the invention as [0.002 S / m, 0.2 S / m, 0.01 S / m, 0°, 0°, 30°], and the parameters are continuously distributed anisotropically along the spherical shell, while the background spherical layer is an isotropic medium, the overlying layer is the crust [0.01 S / m, 0.01 S / m, 0.01 S / m, 0°, 0°, 0°] with a thickness of 80 km, the air layer is [10 -8 S / m, 10 -8 S / m, 10 -8 S / m, 0°, 0°, 0°] with a thickness of 25,484 km, and the core is [1 S / m, 1 S / m, 1 S / m, 0°, 0°, 0°] with a radius of 6071 km. Orthogonal non-planar wave sources (with an external source coefficient of 100 nT) are applied on the outer layer of the mesh, and the right-hand side source term integral is calculated using triangular Gaussian integration. Then, the mass matrix and stiffness matrix of the overall element are constructed according to the relationship between the local numbering and global numbering of the mesh elements. After that, the corresponding large-scale simulation equations are constructed in a frequency loop, the solution method is selected. In this embodiment, the direct solution method is used to solve the equations. After obtaining the edge electric fields of all mesh elements at the current frequency, the magnetotelluric response at the measuring station is calculated through an interpolation function. The loop calculation is repeated until all 10 observation frequencies are calculated, and the magnetotelluric apparent resistivity and phase response of all frequencies and all measuring stations are output as the simulation response for comparative analysis.

[0150] In this embodiment, the simulation response in the three-dimensional spherical coordinate system is compared with the analytical solution of one-dimensional layered magnetotellurics (refer to Pek and Santos, 2002). AsFigure 5 The maximum relative error between the three-dimensional simulation response calculated by the present invention and the one-dimensional analytical solution within the range of 0.1 Hz to 0.0001 Hz is less than 5%, indicating that the magnetotelluric simulation method of the arbitrary anisotropic conductivity tensor in the spherical coordinate system of the present invention is effective. At the same time, there is an obvious gap between the three-dimensional spherical coordinate response above 0.0001 Hz and the one-dimensional analytical solution. This is because the magnetotelluric response above 0.0001 Hz is severely affected by the earth's curvature, and the analytical solution based on the plane wave assumption can no longer meet the simulation requirements, indicating that the magnetotelluric simulation method of the arbitrary anisotropic conductivity tensor in the spherical coordinate system of the present invention is necessary and effective for the implementation of the long-period magnetotelluric observation array at the intercontinental scale.

[0151] The above-described specific embodiments have detailed the technical solutions and beneficial effects of the present invention. It should be understood that the above is only the most preferred embodiment of the present invention and is not used to limit the present invention. Any modifications, supplements, equivalent replacements, etc. made within the scope of the principles of the present invention shall be included in the protection scope of the present invention.

Claims

1. A magnetotelluric simulation method for an arbitrary anisotropic conductivity tensor in spherical coordinates, characterized in that, It includes the following steps: Determine the coordinates and observation frequencies of the intercontinental-scale measurement point array according to the magnetotelluric simulation target; Construct an earth background model and discretize it using unstructured tetrahedral meshes to obtain unstructured meshes; Assign electrical properties to the unstructured meshes using a spherical coordinate arbitrary anisotropic conductivity tensor, including: The expression in spherical coordinates is given as: ； where the superscript s indicates that the physical quantity is in spherical coordinates, i.e., represents the anisotropic conductivity tensor in spherical coordinates. The right side of the formula is in expanded form, where the subscript of each component indicates the scalar conductivity corresponding to it along different directions of the spherical coordinate system. The tensor will be further transformed into three principal conductivity in spherical coordinates and three principal rotation angles in spherical coordinates for description: ； where the superscript T represents the transpose of the matrix, and , , are the three principal axes in spherical coordinates , and corresponding spatial Euler rotation matrices, expressed as: ， ， ； The Ohm's law is satisfied by the conductivity tensor with arbitrary anisotropy in spherical coordinates , where the superscript s represents the physical quantity in spherical coordinates, represents the current density in spherical coordinates, represents the electric field in spherical coordinates; The equivalent form in Cartesian coordinates of the conductivity tensor with arbitrary anisotropy in spherical coordinates is given , where the superscript c represents the physical quantity in Cartesian coordinates, represents the current density in Cartesian coordinates, represents the electric field in Cartesian coordinates, and the equivalent form in Cartesian coordinates of the conductivity tensor with arbitrary anisotropy in spherical coordinates is: ； wherein represents the anisotropic conductivity tensor in the Cartesian coordinate system, and are rotation matrices in two directions, and their expressions are as follows: ， ； Based on this, the electrical properties of the assigned unstructured grid include ; Apply an orthogonal non-planar wave field source to the outer layer of the unstructured grid and calculate the source term integral of the non-planar wave field source, including: specifying two types of orthogonal magnetic fields H As the field source of the external boundary, assuming air insulation, the magnetic field is represented by the negative gradient of the magnetic scalar potential U as: ； Superscript ext Indicates the physical quantity on the outer boundary of the grid, i.e., Indicates the magnetic field on the outer boundary, Indicates the magnetic scalar potential on the outer boundary, Indicates taking the gradient of the physical quantity, Is the magnetic permeability in vacuum, and the magnetic scalar potential in spherical coordinates U Is expressed as spherical harmonics: ； wherein represents the imaginary unit, is the radius of the Earth, represents the external source coefficient, exp represents the exponential function calculation, P is the associated Legendre polynomial, m and n are the degree and order of the polynomial respectively. In the present invention, the field distribution in the quiet Sun state is taken, i.e., m = 0, n = 1. When r is taken as 5 times the radius of the Earth, the first type of field source at the outer boundary is obtained as: ； Among them, and as well as are the magnetic fields in three directions in the spherical coordinate system. According to the conversion relationship between spherical coordinates and Cartesian coordinates, the first type of field source in the Cartesian coordinate system ( , , ) is: ； The second type of field source ( , , ) takes the orthogonal field of the first type of field source as: ； Among them, and as well as are the magnetic fields in three directions of the first type of field source in the Cartesian coordinate system, and as well as are the magnetic fields in three directions of the second type of field source in the Cartesian coordinate system; Use Gaussian integration to calculate the integral of the right-hand source term. The function inside the integral term is regarded as a function related only to the spatial position and is expressed as F. The source term integral is obtained using the integration nodes of the spatial triangle; ； where the superscripts 1, 2 represent the types of field sources, i.e., and represent the source term integrals corresponding to the first type of field source and the second type of field source respectively, A is the weighting coefficient, and the subscript k represents the integration node index; Calculate the mass matrix and stiffness matrix of the elements in the unstructured meshes, and form the overall element matrix according to the numbering relationship; Construct the corresponding simulation equations according to the overall element matrix and the source term integral and at the current observation frequency; Solve the simulation equations to obtain the electric fields of all element edges, and calculate the magnetotelluric response data based on the electric fields through interpolation functions to achieve magnetotelluric simulation.

2. The magnetotelluric simulation method for arbitrarily anisotropic conductivity tensors in spherical coordinate system according to claim 1, wherein Calculate the mass matrix and stiffness matrix of the elements in the unstructured meshes, and form the overall element matrix according to the numbering relationship, including: The electric field double curl equation satisfied by the geomagnetic field is analyzed using the vector finite element method to obtain the mass matrix of element ele and the stiffness matrix as follows: ； ； Among them, the subscript i,j represents the local number of the edge in the element ele, v represents the volume of the element ele, represents the j th edge corresponding to the i th vector interpolation basis function, represents the j th vector interpolation basis function of the edge, × represents taking the curl of the vector interpolation basis function, represents the arbitrary anisotropic conductivity tensor in the spherical coordinate system. Each edge of the element has a global number, so the overall element matrices S and M are formed.

3. The magnetotelluric simulation method for arbitrarily anisotropic conductivity tensors in spherical coordinate system according to claim 1, characterized in that Construct the corresponding simulation equations according to the overall element matrix and the source term integral and at the current observation frequency, including: According to the current observation frequency f Calculate the angular frequency , is the pi, and the corresponding complex frequency coefficient is , where represents the imaginary unit. Let the matrix , then the simulation equations under the current frequency are calculated as follows: ； Where S is the overall stiffness matrix in the overall element matrix, M is the overall mass matrix, the elements in the vector b in the right-hand term are the source term integrals, and only the corresponding positions of the outer boundary edges have values in the vector b, and the rest are 0. The superscripts 1 and 2 represent the types of field sources, and E represents the electric field.

4. The magnetotelluric simulation method for arbitrarily anisotropic conductivity tensors in spherical coordinate systems according to claim 3, characterized in that, Solve the simulation equations to obtain the electric fields of all element edges, including: When the number of simulation units is small and the computing resources are sufficient, a direct solver is used for solution. Its principle is to perform LU decomposition on K of the simulation equations, and back substitution for different right-hand sides of the equations to obtain the electric fields E of all edges 1 and E 2 .

5. The magnetotelluric simulation method for arbitrarily anisotropic conductivity tensors in a spherical coordinate system according to claim 3, characterized in that Solve the simulation equations to obtain the electric fields of all element edges, including: When the number of simulation elements is large and the computing resources are limited, use the iterative method for parallel solution. Use the flexible generalized minimum residual iterative method based on the auxiliary space preconditioner and the block diagonal right preconditioner. First, transform the complex equations presented in the original simulation equations into an equivalent real form, remove the superscripts of different field sources, and give the general form: ； where the superscript re represents the real part of the physical quantity, and the superscript im represents the imaginary part of the physical quantity. Denote the equation in general form as K r ·E r =b r , and the block diagonalization matrix P is: ； Multiply the inverse matrix of P on the right side of the coefficient matrix of the general form equation, then the original equation is K r P -1 ·u = b r , where the vector u is an intermediate vector. The solution process is divided into inner and outer loops. The inner loop uses the auxiliary space preconditioner and the conjugate gradient method to solve the vector , q is the iterative update step vector of the intermediate vector u. The outer loop uses the block diagonal matrix P corresponding to the block diagonal right preconditioner and the generalized minimum residual method to solve the electric field E r .

6. The magnetotelluric simulation method for an arbitrarily anisotropic conductivity tensor in a spherical coordinate system according to claim 1, wherein Calculate the magnetotelluric response data based on the electric fields through interpolation functions to achieve magnetotelluric simulation, including: Calculating the magnetotelluric response data based on the electric fields through interpolation functions includes the magnetotelluric tensor impedance response, apparent resistivity, and phase at the current observation frequency.

Images