Any anisotropic conductivity tensor magnetotelluric simulation method of spherical coordinate system

Abstract

The invention discloses an arbitrary anisotropic conductivity tensor magnetotelluric simulation method of a spherical coordinate system, which belongs to the technical field of geophysical forward modeling, and comprises the following steps: determining the coordinate and observation frequency of a measuring point array; constructing an earth background model and discretizing by adopting a non-structural tetrahedral mesh to obtain a non-structural mesh; assigning a value to the electrical property of the unstructured grid by using any anisotropic conductivity tensor of the spherical coordinates; applying an orthogonal non-planar wave field source to an outer grid of the unstructured grid, and calculating a source item integral of the non-planar wave field source; calculating a mass matrix and a stiffness matrix of the units in the unstructured grid, and forming an overall unit matrix according to the numbering relationship; constructing a corresponding simulation equation set according to the overall unit matrix, the source item integral and the current observation frequency; and solving the simulation equation set to realize magnetotelluric simulation so as to solve the technical problem that the existing magnetotelluric simulation technology cannot flexibly simulate the magnetotelluric response of any anisotropic medium at intercontinental scale.

Description

Technical Field

[0001] The invention belongs to the technical field of geophysical forward modeling, and in particular relates to a spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation method. Background Art

[0002] The activities and evolution of the earth not only bring resources to human beings, but also bring natural disasters. The study of scientific issues related to the earth is of great significance to the sustainable development of mankind. Geophysics is an effective method to detect the distribution of the earth's material structure through physical fields. It has been widely used in disaster monitoring, resource exploration, and research on earth science issues. The detection methods for the mantle scale in the field of geophysics mainly include seismic observation and magnetotelluric sounding. Among them, magnetotelluric sounding uses natural field sources for ground detection. It has many advantages such as large exploration depth and low cost. It is a sharp weapon for deep earth detection. Nowadays, in order to give full play to the detection advantages of magnetotelluric method and realize accurate imaging of the electrical structure of the deep earth, a number of intercontinental-scale magnetotelluric sounding programs have been deployed internationally. Intercontinental-scale magnetotelluric observations can realize the overall imaging of underground media on a large scale and reduce multi-solutions, but it also brings new challenges to the inversion of data.

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

[0004] Unstructured tetrahedron can flexibly fit structures of various scales and is the preferred method for global magnetotelluric three-dimensional forward modeling. Secondly, with the increase of simulation scale, the form of field source has changed significantly. Traditional magnetotelluric field sources are assumed to come from the ionosphere system far away from the earth's surface (4-10 times the earth's radius). For small-scale survey areas, it can be considered that uniform plane polarized waves are vertically incident on the earth's medium. However, when the survey area is intercontinental, this assumption will no longer apply. At this time, a spherical coordinate system field source that conforms to the ionosphere structure should be established. Based on the finite difference algorithm, the magnetotelluric simulation of isotropic medium in spherical coordinate system is realized, and the response differences between plane wave field sources and non-uniform electric field sources are compared. It is believed that the field source will affect the long-period data. Finally, the inhomogeneity of the mantle leads to electrical anisotropy, which will significantly affect the magnetotelluric observation signals. The existing magnetotelluric methods in spherical coordinates are based on the isotropy assumption and easily lead to erroneous results. There are few studies on simulating electrically anisotropic media in spherical coordinates. One method is to use the conductivity anisotropy tensor in the Cartesian coordinate system to simulate the geomagnetic depth sounding response in the spherical coordinate system, but this obviously has serious deviations, and its anisotropy tensor elements 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 data processing needs of intercontinental magnetotelluric arrays in deep earth observations, it is urgent to develop a three-dimensional magnetotelluric simulation method for arbitrary anisotropic media in the spherical coordinate system to serve the future large-scale array data inversion and deep earth science research. Summary of the invention

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

[0007] To achieve the above-mentioned purpose of the invention, an embodiment provides a spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation method, comprising the following steps: Determine the coordinates and observation frequency of the measurement point array according to the magnetotelluric simulation target; Construct the earth background model and use unstructured tetrahedral mesh to discretize the unstructured mesh; Apply spherical coordinate arbitrary anisotropic conductivity tensor to assign electrical properties to unstructured grids; Apply an orthogonal non-plane wave field source to the outer grid of the unstructured grid and calculate the source term integral of the non-plane wave field source; 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; According to the overall unit matrix and source term integration and according to the current observation frequency, the corresponding simulation equation group is constructed; The electric field of all unit edges is obtained by solving the simulation equations, and the magnetotelluric response data is calculated through the interpolation function based on the electric field to realize magnetotelluric simulation.

[0008] Preferably, the electrical properties of the unstructured grid are assigned using an arbitrary anisotropic conductivity tensor in spherical coordinates, including: The expression in spherical coordinates is given as: ; The superscript s Represents the physical quantity in spherical coordinates, that is represents the anisotropic conductivity tensor in spherical coordinates, and the right side of the formula is The expanded form of , where the subscript of each component indicates the corresponding scalar conductivity along different spherical coordinate directions, and the spatial Euler rotation is used to further transform the tensor into the conductivity of the three spherical coordinate principal axes. and the three spherical coordinate principal axis rotation angles To describe: ; The superscript T indicates the transpose of the matrix, and , , are the three principal axes in spherical coordinates , as well as The corresponding spatial Euler rotation matrix is ​​expressed as: , , ; The arbitrary anisotropic conductivity tensor in spherical coordinates satisfies Ohm's law , here superscript s Represents the physical quantity in spherical coordinates, represents the current density in spherical coordinates, represents the electric field in spherical coordinates; The Cartesian equivalent form of the arbitrary anisotropic conductivity tensor in spherical coordinates is given. , here superscript c Represents physical quantities in the Cartesian coordinate system. represents the current density in Cartesian coordinates, Representing the electric field in Cartesian coordinates, the Cartesian equivalent form of the arbitrary anisotropic conductivity tensor in spherical coordinates is obtained as: ; in represents the anisotropic conductivity tensor in Cartesian coordinates, and The rotation matrix in two directions is expressed as follows: , ; Based on this, the electrical properties of the assigned unstructured grid include .

[0009] Preferably, applying an orthogonal non-plane wave field source to the outer grid of the unstructured grid comprises: Given two types of orthogonal magnetic fields H Assuming air insulation as the field source on the external boundary, the magnetic field is expressed as the magnetic scalar potential U The negative gradient of is expressed as: ; Superscript ext Indicates that the physical quantity is on the outer boundary of the grid, that is represents the magnetic field at the outer boundary, represents the magnetic scalar potential at the outer boundary, It means to find the gradient of a physical quantity. is the magnetic permeability in vacuum, the magnetic scalar potential in spherical coordinates U Expressed as spherical harmonics: ; in represents the imaginary unit, is the radius of the Earth, represents the exogenous coefficient, exp represents the exponential function calculation, P To associate the Legendre polynomial, m and n are the degree and order of the polynomial respectively. The present invention takes the field distribution under the quiet day state, that is, m=0, n=1, and takes r as 5 times the radius of the earth to obtain the first type of field source at the outer boundary: ; in, and as well as is the magnetic field 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 ( , , )for: ; The second type of source ( , , ) Take the orthogonal field of the first type of field source as: ; in, and as well as is the magnetic field in three directions of the first type of field source in the Cartesian coordinate system, and as well as is the magnetic field in three directions of the second type field source in the Cartesian coordinate system.

[0010] Preferably, calculating the source term integral of the non-plane wave field source comprises: Gaussian integral is used to calculate the integral of the source term on the right. The integral term is regarded as a function related only to the spatial position and expressed as F. The integral node of the spatial triangle is used to obtain the source term integral: ; Among them, the superscripts 1 and 2 indicate the type of field source, that is, and denote 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.

[0011] Preferably, the mass matrix and stiffness matrix of the cells in the unstructured grid are calculated, and the overall cell matrix is ​​formed according to the numbering relationship, including: The vector finite element method is used to analyze the electric field double curl equation satisfied by the magnetotelluric field, and the mass matrix of the unit ele is obtained. and the stiffness matrix for: ; ; Among them, the subscript i,j represents the local number of the edge in the unit ele, v represents the volume of the unit ele, Indicates j The opposite edge i The vector interpolation basis function corresponding to the edge, Indicates j The vector interpolation basis function of the edge, × means to find the curl of the vector interpolation basis function. Represents an arbitrary anisotropic conductivity tensor in spherical coordinates. Each element edge has a global number, thus forming the overall element matrices S and M.

[0012] Preferably, the corresponding simulation equations are constructed according to the overall unit matrix and the source term integral and according to the current observation frequency, including: According to the current observation frequency f Calculating angular frequency , is the ratio of pi, and the corresponding complex frequency coefficient is ,in Denotes the imaginary unit, let the matrix , then the simulation equations at the current frequency are calculated as follows: ; Among them, S is the overall stiffness matrix in the overall unit matrix, M is the overall mass matrix, the elements in the vector b in the right-hand term are the source term integrals, 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 indicate the type of field source, and E represents the electric field.

[0013] Preferably, solving the simulation equations to obtain the electric fields of all unit edges includes: When the number of simulation units is small and the computing resources are sufficient, a direct solver is used to solve the problem. The principle is to perform LU decomposition on the simulation equation group K and perform back substitution on the right-hand side of different equations to obtain the electric field E of all edges. 1 、E 2 .

[0014] Preferably, solving the simulation equations to obtain the electric fields of all unit edges includes: When the number of simulation units is large and the computing resources are limited, an iterative method is used for parallel solution. A flexible generalized minimum residual iterative method based on the auxiliary space preconditioner and the block diagonal right preconditioner is used. First, the complex equations presented by the original simulation equation group are converted into equivalent real number forms, and the superscripts of different field sources are removed to give the general form: ; The superscript re Indicates the real part of a physical quantity, the superscript im represents the imaginary part of the physical quantity, and the general form of the equation is K r ·E r =b r , the block diagonalization matrix P is: ; Multiply the inverse matrix of P by the coefficient matrix on the right side of the general form equation, and the original equation becomes K r P -1 u=b r , vector u is an intermediate vector, and the solution process is divided into two layers of inner and outer loops. The inner layer uses the auxiliary space preconditioner and 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 generalized minimum residual method to solve the electric field Er .

[0015] Preferably, the magnetotelluric response data is calculated by an interpolation function based on the electric field to realize magnetotelluric simulation, including: The magnetotelluric response data including the magnetotelluric tensor impedance response, apparent resistivity, and phase at the current observation frequency are calculated based on the electric field through the interpolation function.

[0016] Compared with the prior art, the present invention has the following beneficial effects: The present invention adopts an unstructured tetrahedral grid to flexibly depict the real earth environment and finely depict the coastline. At the same time, its hypotenuse characteristics can simulate irregular anomalies (such as subducting plates, folded orogenic belts, mantle plumes, etc.). The present invention adopts spherical harmonic functions to construct non-uniform orthogonal field sources to break the problem that traditional plane wave field sources cannot adapt to the curvature of the earth. At the same time, the present invention proposes a new conductivity tensor of arbitrary anisotropic media in spherical coordinates. The invented tensor is continuously distributed along the spherical coordinates, which solves the problem of discontinuity of the conductivity tensor in the traditional Cartesian coordinates. At the same time, the invention gives an equivalent form of the tensor in Cartesian coordinates based on Euler rotation, so that the unstructured vector finite element method can be applied to construct a group of magnetotelluric simulation equations for arbitrary isotropic media in spherical coordinates. Finally, the invention provides a solution method for the group of electromagnetic simulation equations. According to a specific embodiment provided by the present invention, it can be seen 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 intercontinental-scale magnetotelluric data inversion and deep earth research. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0018] Figure 1 It is a flow chart of a method for simulating arbitrary anisotropic conductivity tensor magnetotelluric in a spherical coordinate system provided by an embodiment; Figure 2 Schematic diagram of the Euler rotation of an arbitrary anisotropic conductivity tensor in spherical coordinates provided in an embodiment; Figure 3 is a schematic diagram of transforming a spherical coordinate system to a Cartesian coordinate system provided by an embodiment; Figure 4 is a schematic cross-sectional view of a model provided in an embodiment and a schematic view of an unstructured tetrahedron mesh (with the air layer removed); Figure 5It is a comparison diagram of the spherical coordinate and analytical solution magnetotelluric responses in the embodiment. DETAILED DESCRIPTION

[0019] To make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific implementation methods described herein are only used to explain the present invention and do not limit the scope of protection of the present invention.

[0020] like Figure 1 As shown, the embodiment provides a spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation method, comprising the following steps: S1, determine the coordinates and observation frequency of the measurement point array according to the magnetotelluric simulation target.

[0021] When actually implementing intercontinental-scale magnetotelluric observations, a measurement network consisting of measurement stations will be established in the measurement area according to the magnetotelluric simulation targets. Each measurement station has corresponding longitude and latitude coordinates, and 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, e is a natural constant of approximately 2.718281828459045), n frequency points (0.1Hz-0.00001Hz) are selected within the period. Therefore, the measurement point array coordinates corresponding to the simulation target are the longitude and latitude coordinates corresponding to the measurement station, and the observation frequency is n frequency points.

[0022] S2, construct the earth background model and use unstructured tetrahedral mesh to discretize to obtain unstructured mesh.

[0023] The shape of the real earth is a complex irregular ellipsoid. Different ellipsoid parameters can be used to construct the earth background model according to the accuracy of the simulation. 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 ocean are portrayed according to the longitude and latitude of the coastline. The earth's interior uses layered spherical shells to describe the crust, mantle, core, etc. If there are specific earth structures that need to be simulated, such as subducting plates and mantle plumes, they can also be simply added to the background grid. The outside of the earth is considered to be an air layer, and a spherical shell of 5 times the radius of the earth is generally used to portray 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 to 6: (1); (2); (3); (4); (5); (6); in x , y , z corresponds to Cartesian coordinates, and r is the radial distance from the sphere to the center of the sphere; is the co-latitude, ranging from 0° to 180°; It is longitude, ranging from 0° to 360°.

[0024] With the Cartesian coordinate model, unstructured tetrahedron grid is used for discretization. The unstructured tetrahedron can flexibly encrypt the coastline, measuring points and places with drastic electrical changes to improve the accuracy of the simulation response.

[0025] S3, use the spherical coordinate arbitrary anisotropic conductivity tensor to assign electrical properties to the unstructured grid.

[0026] The anisotropic conductivity tensor can simultaneously describe the ability of anisotropic and isotropic electrical properties. The present invention provides an expression in spherical coordinates as follows: (7); The superscript s Represents the physical quantity in spherical coordinates, that is represents the anisotropic conductivity tensor in spherical coordinates, and the right side of formula (7) is The expanded form of , where the subscript of each component indicates the corresponding scalar conductivity along different spherical coordinate system directions. From formula (7), it can be seen that the conductivity tensor is continuous along the sphere, and this tensor is positive definite symmetric, satisfying the energy loss law. In order to facilitate simulation, the present invention will further convert the tensor into the conductivity of the three spherical coordinate principal axes by using the spatial Euler rotation. (unit: Siemens per meter, S / m) and the three spherical coordinate main axis rotation angles (unit °), that is: (8); In formula (8), the superscript T represents the transpose of the matrix, and , , are the three principal axes in spherical coordinates , as well as The corresponding spatial Euler rotation matrix is ​​expressed as: , , (9); Figure 2 The Euler rotation method of the spherical coordinate arbitrary anisotropic conductivity tensor of the present invention is described, and the main axis conductivity is rotated around Axis--with rotation Axis--with rotation Axis rotation to new coordinate system , and 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 original spherical coordinate system The equivalent current density under this condition, at this time, the arbitrary anisotropic conductivity tensor in the spherical coordinate system satisfies Ohm's law , here superscript s It also expresses the physical quantity in spherical coordinates, represents the current density in spherical coordinates, represents the electric field in spherical coordinates.

[0027] Since the unstructured grid is established in the Cartesian coordinate system, the present invention provides the Cartesian coordinate equivalent form of any anisotropic conductivity tensor in the spherical coordinate system (satisfying ), here the superscript c Represents physical quantities in the Cartesian coordinate system. represents the current density in Cartesian coordinates, represents the electric field in Cartesian coordinates, such as Figure 3 Describes the transformation of the tensor's spherical coordinate system to the Cartesian coordinate system, that is, the tensor in the spherical coordinate system is first transformed into the Cartesian coordinate system. Axis rotation Angle, so that the rotation The axis is parallel to the z axis, and then revolves around the new Axis rotation Angle, so The axis is parallel to the y-axis. At this time, the orthogonal relationship shows that the axis Parallel to the x-axis, the corresponding electric field is rotated in the same way, and then the current density in the parallel 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 Cartesian equivalent form of the arbitrary anisotropic conductivity tensor in the spherical coordinate system is obtained: (10); In formula (10), represents the anisotropic conductivity tensor in Cartesian coordinates, and The rotation matrix in two directions is expressed as follows: , (11); The angle in the transformation is fixed and unique for each unstructured tetrahedral unit, so it does not need to be given during modeling. At the same time, its equivalent form for isotropic media is consistent with the spherical coordinate form, so this transformation only needs to be applied to the units of anisotropic media during finite element calculations. The electrical properties of the unstructured tetrahedral unit that need to be assigned are .

[0028] S4, applies an orthogonal non-plane wave field source to the outer grid of the unstructured grid and calculates the source term integral of the non-plane wave field source.

[0029] The total field on the outermost boundary of the outer 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 by the present invention is a spherical shell that is 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 magnetotelluric field is unknown, and its response function usually requires an assumed source to describe it. The observation of anisotropic media requires the calculation of the tensor impedance response function, so it is necessary to assume two types of orthogonal sources. The present invention intends to give two types of orthogonal magnetic fields. H Assuming air insulation as the field source on the external boundary, the magnetic field can be expressed as U The negative gradient of is expressed as: (12); In formula (12), the superscript ext Indicates that the physical quantity is on the outer boundary of the grid, that is represents the magnetic field at the outer boundary, represents the magnetic scalar potential at the outer boundary, It means to find the gradient of a physical quantity. is the magnetic permeability in vacuum, the magnetic scalar potential in spherical coordinates U It can be expressed as a spherical harmonic function, namely: (13); in, represents the imaginary unit, is the radius of the Earth, represents the exogenous coefficient, exp represents the exponential function calculation, P To associate the Legendre polynomial, m and n are the degree and order of the polynomial respectively. The present invention takes the field distribution under the quiet day state, that is, m=0, n=1. Taking r as 5 times the radius of the earth, the first type of field source at the outer boundary is: (14); in, and as well as is the magnetic field in three directions in the spherical coordinate system. According to the conversion relationship between spherical coordinates and Cartesian coordinates, we can know that: (15); (16); (17); in, and as well as is the magnetic field in three directions in the Cartesian coordinate system, so the first type of field source in the Cartesian coordinate system ( , , )for: (18); Therefore, the second type of field source ( , , ) Take the orthogonal field of the first type of field source as: (19); When constructing the simulation equations, the source term on the right side is the integral of the field source on the boundary surface. ,in, is the unit vector interpolation basis function, is the normal vector of the face, Indicates that the integral domain is limited to the boundary surface. The present invention uses Gaussian integral to calculate the integral of the right-hand source term. The integral term can be regarded as a function related only to the spatial position and expressed as F. By applying the integral node of the spatial triangle, the source term integral can be obtained as: (20); Among them, the superscripts 1 and 2 indicate the type of field source, that is, and denote 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, with values ​​of 1, 2, 3, and 4, indicating that each triangle has four Gaussian integration nodes.

[0030] S5, calculates the mass matrix and stiffness matrix of the elements in the unstructured grid, and forms the overall element matrix according to the numbering relationship.

[0031] By using the vector finite element method to analyze the electric field double curl equation satisfied by the magnetotelluric field, the mass matrix of the unit ele can be obtained. and the stiffness matrix for: (twenty one); (twenty two); Among them, the subscript i,jrepresents the local number of the edge in the unit ele, v represents the volume of the unit ele, Indicates j The opposite edge i The vector interpolation basis functions corresponding to the edges, each tetrahedron has 6 edges, and the number of edges of each edge is also 6 and includes itself, Indicates j The vector interpolation basis function of the edge, symbol × indicates that the curl of the vector interpolation basis function is calculated. The edges of each unit have a global number, so the overall unit matrices S and M can be formed.

[0032] S6, construct the corresponding simulation equation set according to the overall unit matrix and source term integration and the current observation frequency.

[0033] According to the current observation frequency f Calculating angular frequency , is the ratio of pi, and the corresponding complex frequency coefficient is ,in Denotes the imaginary unit, let the matrix , then the large simulation equations at the current frequency are calculated as follows: (twenty three); Among them, S is the overall stiffness matrix in the overall unit matrix, M is the overall mass matrix, the elements in the vector b in the right-hand term are the source term integrals. From S4, we can see 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 indicate the type of field source, and E represents the electric field.

[0034] S7, solve the simulation equations to obtain the electric field of all unit edges.

[0035] In order to solve the large simulation equation group, namely, equation (23), the present invention proposes to flexibly adopt different solution methods according to the equation scale and computing hardware conditions. When the number of simulation units is small and the computing resources are sufficient, a direct solver is used for solving. The principle is to perform LU decomposition on K of the simulation equation group, and back-substitute the right-hand side terms of different equations to obtain the electric field E of all edges. 1 、E 2 .

[0036] 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 iteration method based on the auxiliary space preconditioner and the block diagonal right preconditioner. To apply this method, the complex equation presented in the original formula (23) needs to be converted into an equivalent real number form. Here, the superscripts of different field sources are removed to give a general form: (twenty four); The superscript re Indicates the real part of a physical quantity, the superscript im represents the imaginary part of the physical quantity, and equation (24) is denoted as K r ·E r =b r , the block diagonalization matrix P is: (25); Multiply the inverse matrix of P by the coefficient matrix on the right side of equation (24), and the original equation becomes K r P -1 u=b r , vector u is an intermediate vector, and the solution process is divided into two layers of inner and outer loops. The inner layer uses the auxiliary space preconditioner and 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 diagonalization matrix P corresponding to the block diagonal right preconditioner and the flexible generalized minimum residual method to solve the electric field E r .

[0037] S8, calculating the magnetotelluric response data through an interpolation function based on the electric field to realize magnetotelluric simulation.

[0038] Solving the large simulation equations, we get the electric field E on the edges of all unstructured tetrahedral units. Applying the interpolation basis function of the unit where the measuring point is located, we can calculate the electric field Ex, Ey, Ez of the measuring point. The corresponding magnetic field can be calculated by Faraday's law of electromagnetic induction. ; Electric field in spherical coordinates E and magnetic field H The electric and magnetic fields can be calculated using Cartesian coordinates: (26); (27); (28); (29); (30); (31); In magnetotelluric observations, the stations are arranged orthogonally in the north-east direction. Therefore, the observation system of the response 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 impedance response Z of the magnetotelluric tensor at the current observation frequency is calculated as: (32); The apparent resistivity is: (33); The phase is: (34).

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

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

[0041] In order to illustrate the technical effect of the above-mentioned spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation method, a specific experimental example is given. Figure 1 The spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation process shown in the figure is implemented. First, it is confirmed that the simulation goal of the embodiment is to verify whether the magnetotelluric response of arbitrary anisotropic media in three-dimensional spherical coordinates is effective and to illustrate the necessity of using spherical coordinates for intercontinental-scale long-period magnetotelluric by comparing the analytical solution of one-dimensional layered anisotropic magnetotelluric response. Therefore, according to the simulation goal, 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 from 0.1 Hz to 0.00001 Hz as the observation frequency, and then as shown in Figure 4 As shown in the figure, the earth background model is constructed (an ideal sphere is selected, in which the radius of the earth is 6371 km). In order to compare with the one-dimensional layered response, the interior of the earth is simply divided into several layers, including an 80 km thick crust and a 220 km thick mantle anisotropic layer. The remaining internal structure is equivalent to a core structure with a radius of 6071 km. Then, it is converted into Cartesian coordinates and discretized using an unstructured tetrahedral grid. The measurement point positions are locally encrypted to obtain a grid of 4685766 unstructured tetrahedral units in total. Then, the electrical properties of the grid are assigned. It is assumed that the mantle has an anisotropic layer with a thickness of 220 km due to its inhomogeneity. The corresponding spherical coordinate anisotropic parameters of this layer are set according to the tensor of the invention patent as [0.002S / m, 0.2 S / m, 0.01 S / m, 0°, 0°, 30°], the anisotropic distribution of the parameters along the spherical shell is continuous, while the background spherical layer is an isotropic medium, the overlying layer is the crust [0.01S / m, 0.01S / m, 0.01S / m, 0°, 0°, 0°], the layer thickness is 80km, and the air layer is [10 -8 S / m, 10 -8 S / m, 10 -8S / m,0°,0°,0°], the layer thickness is 25484km, the core is [1S / m, 1S / m, 1S / m,0°,0°,0°], and the radius is 6071km. An orthogonal non-plane wave source is applied to the outer grid (the external source coefficient is 100nT), and the triangle Gaussian integral is used to calculate the right-hand source term integral, and then the mass matrix and stiffness matrix of the overall unit are constructed according to the relationship between the local number and the global number of the grid unit, and then the corresponding large simulation equation group is constructed according to the frequency cycle, and the solution method is selected. In this embodiment, the direct solution method is used to solve the equation group, and the edge electric field of all grid units at the current frequency is obtained. After that, the magnetotelluric response at the measuring station is calculated by the interpolation function, and the 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 simulation responses for comparative analysis.

[0042] This example compares the three-dimensional spherical coordinate system simulation response with the one-dimensional layered magnetotelluric analytical solution (see Pek and Santos, 2002). Figure 5 It is shown that the maximum relative error of the same-dimensional analytical solution of the three-dimensional simulation response calculated by the present invention between 0.1 Hz and 0.0001 Hz is less than 5%, indicating that the spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation method of the present invention is effective. At the same time, there is an obvious gap between the three-dimensional spherical coordinate response greater than 0.0001 Hz and the same-dimensional analytical solution. This is because the magnetotelluric response greater than 0.0001 Hz is seriously affected by the curvature of the earth, and the analytical solution based on the plane wave assumption can no longer meet the simulation requirements. This shows that the spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation method of the present invention is necessary and effective for the implementation of intercontinental-scale long-period magnetotelluric observation arrays.

[0043] The specific implementation methods described above provide a detailed description of 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 intended to limit the present invention. Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for magnetotelluric simulation of arbitrary anisotropic conductivity tensor in spherical coordinate system, characterized in that: The following steps are involved: Determine the coordinates and observation frequency of the measuring point array according to the magnetotelluric simulation target; Construct the earth background model and use the unstructured tetrahedral grid to discretize the unstructured grid; Apply spherical coordinate arbitrary anisotropic conductivity tensor to assign electrical properties to unstructured grids; Apply an orthogonal non-plane wave field source to the outer grid of the unstructured grid and calculate the source term integral of the non-plane wave field source; 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; According to the overall unit matrix and source term integration and according to the current observation frequency, the corresponding simulation equation group is constructed; The electric field of all unit edges is obtained by solving the simulation equations, and the magnetotelluric response data is calculated through the interpolation function based on the electric field to realize magnetotelluric simulation.

2. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 1, characterized in that: Apply spherical coordinate arbitrary anisotropic conductivity tensors to assign electrical properties to unstructured grids, including: The expression in spherical coordinates is given as: ； The superscript s Represents the physical quantity in spherical coordinates, that is represents the anisotropic conductivity tensor in spherical coordinates, and the right side of the formula is The expanded form of , where the subscript of each component indicates the corresponding scalar conductivity along different spherical coordinate directions, and the spatial Euler rotation is used to further transform the tensor into the conductivity of the three spherical coordinate principal axes. and the three spherical coordinate principal axis rotation angles To describe: ； The superscript T indicates the transpose of the matrix, and , , are the three principal axes in spherical coordinates , as well as The corresponding spatial Euler rotation matrix is ​​expressed as: ， ， ； The arbitrary anisotropic conductivity tensor in spherical coordinates satisfies Ohm's law , here 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 Cartesian equivalent form of the arbitrary anisotropic conductivity tensor in spherical coordinates is given. , here the superscript c Represents physical quantities in the Cartesian coordinate system. represents the current density in Cartesian coordinates, Representing the electric field in Cartesian coordinates, the Cartesian equivalent form of the arbitrary anisotropic conductivity tensor in spherical coordinates is obtained as: ； in represents the anisotropic conductivity tensor in Cartesian coordinates, and The rotation matrix in two directions is expressed as follows: ， ； Based on this, the electrical properties of the assigned unstructured grid include .

3. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 1, characterized in that: Apply orthogonal non-plane wave field sources to the outer grid of an unstructured grid, including: Given two types of orthogonal magnetic fields H Assuming air insulation as the field source on the external boundary, the magnetic field is expressed as the magnetic scalar potential U The negative gradient of is expressed as: ； Superscript ext Indicates that the physical quantity is on the outer boundary of the grid, that is represents the magnetic field at the outer boundary, represents the magnetic scalar potential at the outer boundary, It means to find the gradient of a physical quantity. is the magnetic permeability in vacuum, the magnetic scalar potential in spherical coordinates U Expressed as spherical harmonics: ； in represents the imaginary unit, is the radius of the Earth, represents the exogenous coefficient, exp represents the exponential function calculation, P To associate the Legendre polynomial, m and n are the degree and order of the polynomial respectively. The present invention takes the field distribution under the quiet day state, that is, m=0, n=1, and takes r as 5 times the radius of the earth to obtain the first type of field source at the outer boundary: ； in, and as well as is the magnetic field 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 ( , , )for: ； The second type of source ( , , ) Take the orthogonal field of the first type of field source as: ； in, and as well as is the magnetic field in three directions of the first type of field source in the Cartesian coordinate system, and as well as is the magnetic field in three directions of the second type field source in the Cartesian coordinate system.

4. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 3, characterized in that: Computes source term integrals for non-plane wave field sources, including: Gaussian integral is used to calculate the integral of the source term on the right. The integral term is regarded as a function related only to the spatial position and expressed as F. The integral node of the spatial triangle is used to obtain the source term integral: ； Among them, the superscripts 1 and 2 indicate the type of field source, that is, and denote 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.

5. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 1, characterized in that: Calculate the mass matrix and stiffness matrix of the elements in the unstructured grid, and form the overall element matrix based on the numbering relationship, including: The vector finite element method is used to analyze the electric field double curl equation satisfied by the magnetotelluric field, and the mass matrix of the unit ele is obtained. and the stiffness matrix for: ； ； Among them, the subscript i,j represents the local number of the edge in the unit ele, v represents the volume of the unit ele, Indicates j The opposite edge i The vector interpolation basis function corresponding to the edge, Indicates j The vector interpolation basis function of the edge, × means to find the curl of the vector interpolation basis function. Represents an arbitrary anisotropic conductivity tensor in spherical coordinates. Each element edge has a global number, thus forming the overall element matrices S and M.

6. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 1, characterized in that: According to the overall unit matrix and source term integral and the current observation frequency, the corresponding simulation equations are constructed, including: According to the current observation frequency f Calculating angular frequency , is the ratio of pi, and the corresponding complex frequency coefficient is ,in Denotes the imaginary unit, let the matrix , then the simulation equations at the current frequency are calculated as follows: ； Among them, S is the overall stiffness matrix in the overall unit matrix, M is the overall mass matrix, the elements in the vector b in the right-hand term are the source term integrals, 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 indicate the type of field source, and E represents the electric field.

7. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 6, characterized in that: Solve the simulation equations to obtain the electric field of all element edges, including: When the number of simulation units is small and the computing resources are sufficient, a direct solver is used to solve the problem. The principle is to perform LU decomposition on the simulation equation group K and perform back substitution on the right-hand side of different equations to obtain the electric field E of all edges. 1 、E 2 .

8. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 6, characterized in that: Solve the simulation equations to obtain the electric field of all element edges, including: When the number of simulation units is large and the computing resources are limited, an iterative method is used for parallel solution. A flexible generalized minimum residual iterative method based on the auxiliary space preconditioner and the block diagonal right preconditioner is used. First, the complex equations presented by the original simulation equation group are converted into equivalent real number forms, and the superscripts of different field sources are removed to give the general form: ； The superscript re Indicates the real part of a physical quantity, the superscript im represents the imaginary part of the physical quantity, and the general form of the equation is K r ·E r =b r , the block diagonalization matrix P is: ； Multiply the inverse matrix of P by the coefficient matrix on the right side of the general form equation, and the original equation becomes K r P -1 u=b r , vector u is an intermediate vector, and the solution process is divided into two layers of inner and outer loops. The inner layer uses the auxiliary space preconditioner and 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 generalized minimum residual method to solve the electric field E r .

9. The method for spherical coordinate system arbitrary anisotropic conductivity tensor magnetotelluric simulation according to claim 1, characterized in that: Based on the electric field, the magnetotelluric response data is calculated through the interpolation function to realize magnetotelluric simulation, including: The magnetotelluric response data including the magnetotelluric tensor impedance response, apparent resistivity, and phase at the current observation frequency are calculated based on the electric field through the interpolation function.

Images