阅读说明：本技术 基于地空分解策略的大地电磁正演方法及装置、存储介质 (Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium ) 是由 周峰 任政勇 陈煌 邓居智 张志勇 汤井田 李广 蓝泽鸾 于 2020-06-09 设计创作，主要内容包括：本发明公开了基于地空分解策略的大地电磁正演方法及装置、存储介质。该方法包括：包括：依据数字地形数据(DEM)来构建求解区域三维MT的地电几何模型；将地电几何模型离散为一系列四面体单元；获取基于地空分解策略的A-Φ-Ψ耦合势三维MT正演满足的积分弱解形式并对积分弱分解形式进行基于节点离散的转换得到稀疏线性方程组；加载入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域中四面体单元各个节点的矢量位A、标量位Φ,标量位Ψ；根据矢量位A、标量位Φ,标量位Ψ计算出观测点的参数。本发明不仅对开展大规模三维大地电磁法反演研究具有推动作用,对于促进学科发展,提升电磁法资料解释的学术水平具有重要现实意义与价值。(The invention discloses a magnetotelluric forward modeling method and device based on a geospatial solution strategy and a storage medium. The method comprises the following steps: the method comprises the following steps: constructing a geoelectric geometric model for solving the three-dimensional MT of the area according to digital terrain Data (DEM); discretizing the geoelectric geometric model into a series of tetrahedral units; acquiring an integral weak solution form satisfied by forward modeling of the A-phi-psi coupling potential three-dimensional MT based on a geospatial solution strategy, and performing node-discrete-based conversion on the integral weak solution form to obtain a sparse linear equation set; loading an incident field and placing the incident field in a right-end item of a sparse linear equation set, and solving the sparse linear equation to obtain a vector position A, a scalar position phi and a scalar position psi of each node of a tetrahedral unit in a solution area; and calculating parameters of the observation point according to the vector position A, the scalar position phi and the scalar position psi. The method has a promoting effect on the development of large-scale three-dimensional geoelectromagnetic inversion research, and has important practical significance and value on promoting subject development and improving the academic level of electromagnetic method data interpretation.)

1. A magnetotelluric forward modeling method based on a space-time solution strategy is characterized by comprising the following steps: the method comprises the following steps:

step 1: constructing a geoelectric geometric model for solving the three-dimensional MT of the area;

step 2: discretizing the geoelectric geometric model into a series of non-overlapping tetrahedral cells;

and step 3: acquiring an integral weak solution form satisfied by A-phi-psi coupling potential three-dimensional MT forward modeling based on a geospatial solution strategy, and performing node-based discrete conversion on the integral weak solution form by adopting node basis functions and nodes of tetrahedral units to obtain a sparse linear equation set, wherein A represents a magnetic field vector position, phi represents an electric field scalar position, psi represents a magnetic field scalar position, and the nodes are connection hubs between discrete tetrahedral units;

and 4, step 4: loading an incident field and placing the incident field in a right end item of a sparse linear equation set, and solving the sparse linear equation to obtain a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi of each node of a tetrahedral unit in a solving region, wherein the magnetic field vector position A and the electric field scalar position phi of the node in an air region are 0, and the magnetic field scalar position psi of the node in an underground region is 0;

and 5: and calculating the combination of one or more parameters of the electric field and the magnetic field of the observation point according to the magnetic field vector position A, the electric field scalar position phi and the magnetic field scalar position psi.

2. The method of claim 1, wherein: the integral weak solution form in step 3 is as follows:

in the formula, V represents a test function, E represents an electric field, B represents a magnetic field, omega represents a solution area,denotes the boundary of the solution area, Ω₀Is the air space region, omega₁Is a region of a subterranean space,₀which represents the upper top surface of the air layer,₁represents the lower bottom surface of the subsurface region, which is the interface between the subsurface region and the air region, n represents an out-of-boundary normal vector, dv, ds represent a volume differential element and an area differential element, respectively, v represents a volume, s represents an area,the resistance ratio is represented by the ratio of resistance,

3. The method of claim 1, wherein: the sparse linear equation set in step 3 is as follows:

in the formula (I), the compound is shown in the specification,A^f,A^mare all sub-matrices that are sparse,are all sub-matrices formed by the area integral terms of the geospatial interface,all represent a zero matrix, B₀,B₁,B₂,B₃,B₄All right-hand terms of a sparse linear system of equations, X₀、X₁、X₂、X₃、X₄The magnetic field vector position A, the electric field scalar position phi and the magnetic field scalar position psi of all nodes in a solution region are respectively an unknown matrix in a sparse linear equation set.

4. The method of claim 3, wherein: the parameters in the sparse linear system of equations are represented as follows:

X₃＝(Φ₁,Φ₂,…Φ_n)^T，

wherein i, j is 1 to n; m, k is 1 to n₁，

L_i,L_jrespectively representing node basis functions L corresponding to the nodes i and j in the underground area_k,L_mRespectively represents the node basis functions corresponding to the kth node and the mth node in the air region,representing a sub-matrixThe elements corresponding to the middle node i and the node j,represents a sub-matrix B₀The element corresponding to the middle node i is,

5. The method of claim 1, wherein: loading an incident field in a right-end item of a sparse linear equation set in the step 4, solving the sparse linear equation to obtain a magnetic field vector position A and an electric field scalar position phi of each node of a tetrahedral unit in a solved area, wherein the process of the magnetic field scalar position psi is as follows:

①, constructing an A-phi-psi coupling potential three-dimensional MT dual problem based on a geospatial solution strategy, and solving the dual problem to obtain a dual field W_A-W_φ-W_Ψ(ii) a Calculating a right end item of the sparse linear equation set based on the loaded incident field, and solving the sparse linear equation set to obtain a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi of a node corresponding to the current tetrahedron unit dispersion in the solution region;

wherein, W_AA field of parity, W, representing the magnetic field vector bit A_ΦExpressed as electric field scalar potential phiA field of parity, W_ψA field of parity which is the scalar position Ψ of the magnetic field;

② field vector A, electric field scalar phi, and magnetic field scalar psi at each node, and the field vector A, calculated based on step ①_AElectric field scalar potential phi of the field_ΦThe field scalar psi_ψCalculating unit error distributions respectively corresponding to dual problems of the geoelectromagnetic method problem, and calculating relative error unit indicating factors of the tetrahedral units based on the unit error distributions;

③: respectively judging whether the relative error indicating factor of each tetrahedral unit is greater than a preset maximum error, if so, performing refinement dispersion on the tetrahedral units, and returning to the step 3 until an iteration termination condition is met;

the iteration termination condition is as follows: the iteration times reach the preset maximum iteration times or the total number of the tetrahedral units in the solving area reaches the preset maximum unit number;

and finally obtaining a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi on all the tetrahedral unit nodes.

6. The method of claim 5, wherein: the unit error distributions respectively corresponding to the dual problems of the magnetotelluric problem are as follows:

in the formula, e and w represent the cell error distribution corresponding to the magnetotelluric problem and the dual problem, respectively, n represents the out-of-boundary normal vector, I is an imaginary number cell, ω is an angular frequency, ds represents an area differential element, s represents an area, μ is a magnetic permeability,is the surface of a tetrahedral unit, the symbol [ ·]Representing a field in a plane

the relative error unit indicator is as follows:

In the formula, gamma_eA factor is indicated for the relative error cell of a tetrahedral cell, γ is the weighted error of the tetrahedral cell,

7. The method of claim 1, wherein: in step 4, loading an incident field in a right-end item of the sparse linear equation set, dividing the incident field into two polarization directions, and respectively loading the incident field to obtain a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi of each node corresponding to the two polarization directions, wherein the method further comprises the following steps:

and obtaining one or more combinations of impedance tensor, apparent resistivity and phase response based on the magnetic field vector position A and the electric field scalar position phi of each node corresponding to the two polarization directions.

8. The method of claim 1, wherein: before solving the sparse linear equation, acquiring physical parameters of different areas in the geoelectric geometric model;

wherein the physical property parameters include a resistivity parameter, a dielectric constant parameter, and a permeability parameter of each region.

9. An apparatus based on the method of any one of claims 1-8, characterized in that: the method comprises the following steps:

a model construction module: the geoelectricity geometric model is used for constructing a three-dimensional MT of the solution area;

tetrahedral unit discrete module: for discretizing the geoelectric geometric model into a series of non-overlapping tetrahedral cells;

an operation module: and performing forward calculation on the A-phi-psi coupling potential based on the geospatial solution strategy to finally obtain the target parameters of the observation points.

10. A storage medium having a computer program stored thereon, characterized in that: the computer program when executed by a processor implements the steps of the method of any one of claims 1 to 8.

Technical Field

The invention belongs to the technical field of applied geophysics, and particularly relates to a magnetotelluric forward modeling method and device based on a geospatial solution strategy and a storage medium.

Background

The Magnetotelluric (MT) method is a geophysical exploration method for detecting the electrical structure in the earth by using a natural alternating electromagnetic field source, is widely applied to the research in the fields of mineral resource exploration, oil-gas and natural gas resource exploration, geothermal resource exploration, deep geological structure detection and the like, and has good application effect.

At present, the accuracy of the data interpretation of the magnetotelluric method can be directly influenced by any complex field exploration environment, such as an exploration area which is transited from a gentle area to a complex mountain area; secondly, the large-scale three-dimensional data set is emerging, so that the existing MT data interpretation technology faces the difficult problems of large memory consumption and the like. In addition, when solving a large-scale problem based on the existing three-dimensional MT forward solving formula, a system matrix obtained by dispersing the forward solving formula faces the problem of conditional number difference, and the solving efficiency of a linear equation set is seriously influenced. Aiming at the problems, the data interpretation technology of the three-dimensional MT data set under large-scale and arbitrary complex terrains, which is applied to the exploration of deep mineral resources and the exploration of oil gas and natural gas resources, is to be realized, and the forward technology with high speed and high precision is a precondition. Therefore, a high-precision and efficient three-dimensional geoelectromagnetic forward modeling method is needed.

Disclosure of Invention

The invention aims to provide a magnetotelluric forward modeling method based on a geospatial solution strategy, which realizes three-dimensional MT self-adaptive finite element forward modeling based on a geospatial solution strategy A-phi-psi coupling potential, can realize high-efficiency and high-precision forward modeling calculation, overcomes the problems of poor condition number, low solving efficiency, poor precision and the like of the traditional MT solving system, has the capabilities of high efficiency, high precision and high stability, and provides core power for realizing inversion interpretation of measured data of a large-scale three-dimensional magnetotelluric method.

On one hand, the magnetotelluric forward modeling method based on the geospatial solution strategy provided by the invention comprises the following steps:

step 1: constructing a geoelectric geometric model for solving the three-dimensional MT of the area according to digital terrain Data (DEM);

step 2: discretizing the geoelectric geometric model into a series of non-overlapping tetrahedral cells;

and step 3: acquiring an integral weak solution form satisfied by A-phi-psi coupling potential three-dimensional MT forward modeling based on a geospatial solution strategy, and performing node-based discrete conversion on the integral weak solution form by adopting node basis functions and nodes of tetrahedral units to obtain a sparse linear equation set, wherein A represents a magnetic field vector position, phi represents an electric field scalar position, psi represents a magnetic field scalar position, and the nodes are connection hubs between discrete tetrahedral units;

and 4, step 4: loading an incident field and arranging the incident field in a right-end item of a sparse linear equation set, and solving the sparse linear equation to obtain a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi of each node of a tetrahedral unit in a solution area; wherein, the magnetic field vector position A and the electric field scalar position phi of the node in the air region are 0, and the magnetic field scalar position psi of the node in the underground region is 0;

and 5: and calculating the combination of one or more parameters of the electric field and the magnetic field of the observation point according to the magnetic field vector position A, the electric field scalar position phi and the magnetic field scalar position psi.

Further preferably, the form of the integral weak solution in step 3 is as follows:

in the formula, V represents a test function, E represents an electric field, B represents a magnetic field, omega represents a solution area,denotes the boundary of the solution area, Ω₀Is the air space region, omega₁Is a region of a subterranean space,₀which represents the upper top surface of the air layer,₁represents the lower bottom surface of the subsurface region, which is the interface between the subsurface region and the air region, n represents an out-of-boundary normal vector, dv, ds represent a volume differential element and an area differential element, respectively, v represents a volume, s represents an area,the resistance ratio is represented by the ratio of resistance,

denotes the admittance ratio, and μ is the permeability.

Further preferably, the sparse linear equation set in step 3 is as follows:

in the formula (I), the compound is shown in the specification,

A^f,A^mare all sub-matrices that are sparse,are all sub-matrices formed by the area integral terms of the geospatial interface,

all represent a zero matrix, B₀,B₁,B₂,B₃,B₄All right-hand terms of a sparse linear system of equations, X₀、X₁、X₂、X₃、X₄The magnetic field vector position A, the electric field scalar position phi and the magnetic field scalar position psi of all nodes in a solution region are respectively an unknown matrix in a sparse linear equation set.

Further preferably, the parameters in the sparse linear system of equations are expressed as follows:

wherein i, j is 1 to n; m, k is 1 to n₁，

Representing three directions of a cartesian coordinate system,

respectively, the derivation of x is shown,respectively, the derivation of y is shown,

respectively representing the derivation of z [ ·]_xWhich represents the component in the x-direction,

a magnetic field vector A representing the first, second and nth nodes of the underground region in the x direction_x；

A magnetic field vector bit A representing the y direction of the first node, the second node and the nth node of the underground region_y；

A magnetic field vector A representing the first, second and nth nodes of the underground region in the z direction_zMagnetic field vector bit A in x, y, z directions_x、A_y、A_zForming a magnetic field vector bit A; phi₁、Φ₂、Φ_nAn electric field scalar potential representing a first node, a second node, and an nth node of the subsurface region; Ψ₁、Ψ₂、

Representing the first node, the second node, the n-th node of the air zone₁A magnetic field scalar bit of the node, T represents a matrix transposition, n is the total number of nodes of the underground area, n₁The total number of nodes in the air area is defined, and the nodes correspond to the vertexes of the tetrahedral unit;

L_i,L_jrespectively represents the node basis functions corresponding to the nodes i and j in the underground area, respectively represents the node basis functions corresponding to the kth node and the mth node in the air area,representing a sub-matrix

The elements corresponding to the middle node i and the node j,represents a sub-matrix B₀The definition of the element corresponding to the middle node i is similar to that of the other sub-matrices, and is not repeated herein.

The resistance ratio is represented by the ratio of resistance,denotes the admittance ratio, Ω₀Is the air space region, omega₁Is a region of a subterranean space,₀which represents the upper top surface of the air layer,₁the lower bottom surface of the underground area is represented as a boundary surface between the underground area and the air area, n represents an out-of-boundary normal vector, dv and ds represent a volume differential element and an area differential element, respectively, v represents a volume, s represents an area, and μ represents a permeability.

Respectively representing three components, E, of a background vector bit A generated after a polarized field source is excited under a uniform semi-space or laminar dielectric earth model^b,B^bRespectively representing the background electric field and the magnetic field generated after the polarized field source is excited under a uniform semi-space or laminar dielectric electric model. Therefore, the right-end term of the coefficient linear equation system can be calculated by using the formula when the incident field is loaded.

Further preferably, loading an incident field in the right-end item of the sparse linear equation set in step 4, and solving the sparse linear equation to obtain a magnetic field vector position a and an electric field scalar position Φ of each node of the tetrahedral unit in the solution area, wherein the process of the magnetic field scalar position Ψ is as follows:

① construction of A-phi-psi coupling potential three-dimensional MT dual based on geospatial solution strategySolving the problem to obtain a dual field W_A-W_φ-W_Ψ(ii) a Calculating a right end item of the sparse linear equation set based on the loaded incident field, and solving the sparse linear equation set to obtain a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi of a node corresponding to the current tetrahedron unit dispersion in the solution region;

wherein, the dual problem is converted into a linear equation set, the right end item of the linear equation set is obtained by solving a three-component unit source integral expression, and then the linear equation set is solved to obtain a dual field W_A-W_Φ-W_ψ，W_AA field of parity, W, representing the magnetic field vector bit A_ΦA field of parity, W, expressed as the electric field scalar potential phi_ψA field of parity which is the scalar position Ψ of the magnetic field;

② field vector A, electric field scalar phi, and magnetic field scalar psi at each node, and the field vector A, calculated based on step ①_AElectric field scalar potential phi of the field_ΦThe field scalar psi_ψCalculating unit error distributions respectively corresponding to dual problems of the geoelectromagnetic method problem, and calculating relative error unit indicating factors of the tetrahedral units based on the unit error distributions;

③: respectively judging whether the relative error indicating factor of each tetrahedral unit is greater than a preset maximum error, if so, performing refinement dispersion on the tetrahedral units, and returning to the step 3 until an iteration termination condition is met;

the iteration termination condition is as follows: the iteration times reach the preset maximum iteration times or the total number of the tetrahedral units in the solving area reaches the preset maximum unit number;

and finally obtaining a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi on all the tetrahedral unit nodes.

The preset maximum error, the preset maximum iteration number, and the preset maximum unit number may be empirical values obtained by an actual demand experiment.

More preferably, the unit error distributions respectively corresponding to the dual problems of the magnetotelluric problem are as follows:

in the formula, e and w represent the cell error distribution corresponding to the magnetotelluric problem and the dual problem, respectively, n represents the out-of-boundary normal vector, I is an imaginary number cell, ω is an angular frequency, ds represents an area differential element, s represents an area, μ is a magnetic permeability,

is a tetrahedral surface, the symbol [ ·]Representing a field in a plane

A jump value of;

the relative error unit indicator is as follows:

wherein γ ═ ew-

In the formula, gamma_eA factor is indicated for the relative error cell of a tetrahedral cell, γ is the weighted error of the tetrahedral cell,

the maximum value of the weighting errors in all tetrahedral cells.

Further preferably, the loading of the incident field in step 4 is performed in two polarization directions, and the right-hand item of the sparse linear equation set is divided into a magnetic field vector position a, an electric field scalar position Φ, and a magnetic field scalar position Ψ of each node corresponding to the two polarization directions are obtained by loading the incident field in the two polarization directions, and the method further includes:

and obtaining one or more combinations of impedance tensor, apparent resistivity and phase response based on the magnetic field vector position A and the electric field scalar position phi of each node corresponding to the two polarization directions.

And obtaining a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi of each node in the first polarization direction according to the first-third step. The second polarization direction can also be obtained according to the third step (i) and the fourth step (ii) after the magnetic field vector position A, the electric field scalar position phi and the magnetic field scalar position psi of each node. In addition, the second polarization direction can also be directly calculated by the forward modeling grid after the adaptive encryption of the first polarization direction, which takes into account that the difference of the forward modeling grids in different polarization directions in the same solution area is not large.

Further preferably, before solving the sparse linear equation, the method further comprises the step of acquiring physical parameters of different regions in the geoelectrical geometric model;

wherein the physical property parameters include a resistivity parameter, a dielectric constant parameter, and a permeability parameter of each region.

In another aspect, the present invention provides an apparatus of the above method, including:

a model construction module: the geoelectricity geometric model is used for constructing a three-dimensional MT of the solution area;

tetrahedral unit discrete module: for discretizing the geoelectric geometric model into a series of non-overlapping tetrahedral cells;

an operation module: and performing forward calculation on the A-phi-psi coupling potential based on the geospatial solution strategy to finally obtain the target parameters of the observation points.

Furthermore, the present invention provides a storage medium having stored thereon a computer program which, when being executed by a processor, carries out the steps of the above-mentioned method.

Advantageous effects

1. The magnetotelluric forward modeling method based on the geospatial solution strategy is an A-phi-psi coupling potential three-dimensional MT finite element algorithm based on the geospatial solution strategy, the algorithm is carried out in a node finite element space, and compared with a vector finite element of a traditional electric field equation, the phenomenon of an empty solution does not exist. Meanwhile, compared with a vector finite element, the matrix condition number formed by the node finite element is more balanced in value, and the equation system is more stable to solve. The A-phi-psi coupling potential three-dimensional MT finite element algorithm based on the geospatial solution strategy can carry out iterative solution, and vector finite elements of an electric field equation are not suitable for iterative calculation and only suitable for a direct solver. The iterative solver and the direct solver have the most obvious advantages that the memory requirement of a computer is low, and meanwhile, a system matrix is not required to be used for carrying out inversion, so that the time is saved, the forward solving efficiency is improved, and the memory consumption is saved. Therefore, the A-phi-psi coupling potential forward modeling solving system based on the geospatial decomposition strategy provided by the invention can be suitable for fast forward modeling of a large-scale three-dimensional magnetotelluric method, has a promoting effect on inversion research of the three-dimensional magnetotelluric method, and also has important practical significance and value for promoting subject development and improving the academic level of electromagnetic method data interpretation.

2. According to the invention, the numerical precision of the three-dimensional MT finite element mainly depends on mesh subdivision, the traditional technology can meet the design requirement for a simple geoelectric model, for the geoelectric model with complex undulating terrain, the design of the mesh consumes a large amount of time, and an ideal effect is often not obtained.

3. The method comprises the steps of constructing a dual problem of an A-phi-psi coupling potential based on a geospatial solution strategy, and realizing posterior error estimation according to unit weighted errors, thereby realizing an A-phi-psi coupling potential adaptive MT (MT) forward algorithm based on the geospatial solution strategy.

Drawings

FIG. 1 is a schematic diagram of an MT survey. Wherein the content of the first and second substances,₀which represents the upper top surface of the air layer,₁representing the lower floor of the subsurface region, which is undulating, with n representing the outer normal, omega₀Expressed as air space, omega₁Represented as a subterranean space.

FIG. 2 is a forward flow chart of an A-phi-psi coupling potential adaptive three-dimensional magnetotelluric method based on a geospatial solution strategy;

FIG. 3 is a schematic view of a geometric electrical model, wherein (a) is a schematic view of an x-y plane and (b) is a schematic view of an x-z section;

FIG. 4 is a schematic diagram of a three-dimensional terrain model mesh;

FIG. 5 is a graph of ρ from a terrain model at a frequency of 2Hz^xxAnd ρ^xyApparent resistivity contour map, wherein (a) is rho^xxApparent resistivity contour map, and (b) is rho^xyApparent resistivity contour map;

FIG. 6 is a graph of ρ from a terrain model at a frequency of 2Hz^yxAnd ρ^yyApparent resistivity contour map, wherein (a) is rho^yxApparent resistivity contour map, and (b) is rho^yyApparent resistivity contour plot.

Detailed Description

The present invention will be further described with reference to the following examples.

The invention provides a magnetotelluric forward modeling method based on a geospatial solution strategy, aiming at the difficult problems of poor condition number, low solution efficiency, poor precision and the like of the traditional MT solution system. The technology realizes a 3D MT self-adaptive finite element forward modeling technology based on a geospatial solution strategy A-phi-psi coupling potential, can realize high-efficiency and high-precision forward modeling calculation on any three-dimensional complex terrain MT problem, and provides core power for realizing inversion interpretation of measured data of a large-scale three-dimensional magnetotelluric method. It is carried out according to the following principle:

1. regarding the 3D isotropic MT edge value problem:

based on Maxwell equation satisfied by the magnetotelluric method, a vector Helmholtz total field equation for deriving the coupling between the magnetic field vector position A and the electric field scalar position phi is as follows:

wherein the content of the first and second substances,

the admittance ratio is expressed as a function of,representing the impedance ratio, i is an imaginary unit, a represents a magnetic field vector position, Φ and Ψ are an electric field scalar position and a magnetic field scalar position, respectively, σ is the conductivity (reciprocal of the resistivity), ω ═ 2 π f is the angular frequency, f is the frequency, and ≈ 8.854187817 × 10^-12F/m, dielectric constant, μ₀And (4) vacuum magnetic permeability.

As is well known, in the MT band, it can be assumed that there is no displacement current, outside the earth (i.e. air space), the magnetic field B is expressed as:

as can be seen from equation (2), the magnetic field scalar position Ψ satisfies the following equation,

the scalar field position Ψ is solved only in the air region, and the vector field position a and the scalar field position Ψ satisfy the following equation at the air interface according to the field continuity condition,

meanwhile, in the air, the magnetic field vector is located at A and the electric field

The normal direction is not continuous,(J is the current density) and the magnetic field vector bit a vanishes at the ground-space boundary, and therefore,

wherein a ground-air interface is represented, as shown in fig. 1.

The 3D MT problem based on the A-phi-psi of the ground-space region decomposition satisfies the following differential equation system:

in order to ensure the uniqueness of the differential equation set solution, it needs to incorporate the following boundary conditions,

wherein, the boundary between the underground region and the air region is the boundary, and n represents the normal vector outside the boundary.

2. And (3) integrating weak solution forms corresponding to the 3D isotropic MT edge value problem.

The residual vector R can be obtained by the formula (6)_Ψ,R_A,R_ΦThe specific expression of (A) is as follows:

according to the weighted residual method, selecting a test function V, and enabling the weight of the residual error in a calculation region omega (a solving region) to be 0:

substituting the residual term (8) into (9), the differential equation satisfied by the three-dimensional MT of A-phi-psi, which can be decomposed based on the ground-space region, is transformed into the form of integral weak solution, as follows:

in the formula, V represents a test function, E represents an electric field, B represents a magnetic field, omega represents a solution area,

the boundary of the solution area is represented,₀which represents the upper top surface of the air layer,₁representing the lower floor of the subsurface region, dv, ds represent the volume differential element and the area differential element, respectively (v represents the volume and s represents the area).

And converting the integral weak solution form to obtain a sparse linear equation set shown as follows:

in the formula (I), the compound is shown in the specification,A^f,A^mare all sub-matrices that are sparse,are all sub-matrices formed by the area integral terms of the geospatial interface,

all represent a zero matrix, B₀,B₁,B₂,B₃,B₄Are all right-hand terms of a sparse linear system of equations.

X₃＝(Φ₁,Φ₂,…Φ_n)^T(17)

Wherein i, j is 1 to n; m, k is 1 to n₁，

Representing three directions of a cartesian coordinate system,

respectively, the derivation of x is shown,respectively, the derivation of y is shown,respectively, the derivation of z is shown,

a magnetic field vector A representing the first, second and nth nodes of the underground region in the x direction_x；A magnetic field vector bit A representing the y direction of the first node, the second node and the nth node of the underground region_y；

A magnetic field vector A representing the first, second and nth nodes of the underground region in the z direction_zMagnetic field vector bit A in x, y, z directions_x、A_y、A_zForming a magnetic field vector bit A; phi₁、Φ₂、Φ_nAn electric field scalar potential representing a first node, a second node, and an nth node of the subsurface region; Ψ₁、Ψ₂、Representing the first node, the second node, the n-th node of the air zone₁A magnetic field scalar bit of the node, T represents a matrix transposition, n is the total number of nodes of the underground area, n₁The total number of nodes in the air area is defined, and the nodes correspond to the vertexes of the tetrahedral unit;

L_i,L_j,L_k,L_mrepresenting the node basis function, E the electric field, B the magnetic field,the resistance ratio is represented by the ratio of resistance,

denotes the admittance ratio, Ω₀Is the air space region, omega₁Is a region of a subterranean space,₀which represents the upper top surface of the air layer,₁the lower bottom surface of the underground area is represented as a boundary surface between the underground area and the air area, n represents an out-of-boundary normal vector, dv and ds represent a volume differential element and an area differential element, respectively, v represents a volume, s represents an area, and μ represents a permeability. The unknowns of the hybrid solution system are 4n + n₁N is the total number of nodes in the underground region, n₁The total number of nodes in the air zone. In order to solve the integral sparseness of each submatrix, the test function needs to be replaced by a node function of 1 st order, i.e., V is L, and L is a node basis function.

For any non-node position in a tetrahedral unit, calculating a vector position A, a scalar position phi and a scalar position psi as follows:

wherein L is_jAnd representing the node basis function corresponding to the jth node of the tetrahedral unit. If we push to n node sums, then equations (35) - (37) exist: if the node is not located on the tetrahedral unit, then L is at this time_jIs 0.

And finally, calculating the equation set (13) by adopting a GMRES iterative solver to obtain a vector position A and scalar positions phi and psi of each node. It should be noted that the calculated vector position a and scalar position Φ are for a subsurface region, the vector position a and scalar position Φ for a node in an air region are 0, the calculated scalar position Ψ is for an air region, and the scalar position Ψ for a node in a subsurface region is 0. Thus, the vector position A and the scalar positions Φ and Ψ for all nodes in the subsurface region and the air region can be known from the above calculations.

3. Mesh adaptive encryption

In order to realize the target-oriented adaptive finite element technology, an A-phi-psi coupling potential 3D MT dual problem based on a geospatial solution strategy needs to be constructed. Firstly, a simple 'auxiliary artificial source' is introduced into a measuring point area, the 'auxiliary artificial source' is not an actual physical problem and is a mathematical problem constructed artificially, and the aim of the 'auxiliary artificial source' is to encrypt the measuring point area. By using omega_jThe area representing the observation point passes through omega_jAnd (3) constructing a unified dual problem corresponding to the magnetotelluric problem by using a unit source I (1,1,1) of three components of region injection:

whereinW_AA field of parity, W, representing the magnetic field vector bit A_ΦA field of parity, W, expressed as the electric field scalar potential phi_ΨIs a field of parity of the scalar position psi in air,_jthe total field equation and the even field W of the A-phi-psi coupling potential are satisfied by the geoelectromagnetic method by observing the structural form of an equation (38) which is a delta function and t is the number of unit sources of three components_A-W_φ-W_ΨHave the same structure and thus the same finite element linear system matrix. Thus, equation (38) corresponds to a dual problem linear system of equations as:

KY＝C (39)

where K is the stiffness matrix formed by the dual problem, which is consistent with the sparse matrix of equation (13), Y is the dual field at the cell node,

and solving the expression of the three-component unit source integral.

e is the a posteriori error estimate over the tetrahedral unit,is the surface of a tetrahedral unit, the symbol [ ·]Representing a field in a planeThe jump value of (2). The unit error distribution corresponding to the magnetotelluric problem and the dual problem respectively, e and w:

and (3) adopting a weighting idea, wherein the finally obtained weighting error gamma of the tetrahedral unit is as follows:

γ＝|ew| (41)

finally, the relative error cell indicator for the current cell can be estimated

In the formula (I), the compound is shown in the specification,representing the maximum value of the weighting error in all tetrahedral units.

4. And calculating the field value at the observation point.

In order to calculate the actual application of the 3D MT, the electric field and the magnetic field can be calculated by using theoretical formulas in the field, such as an electric field formula:in addition, in order to obtain parameters such as apparent resistivity, phase and the like, the invention needs to load field sources in different polarization directions to obtain the apparent resistivity and the phase in different polarization directions, and the expression for specifically calculating the apparent resistivity and the phase at the observation point is as follows:

in the formula, Z_xx、Z_xy、Z_yx、Z_yyRespectively representing the four components of the impedance tensor Z,respectively representing the electric field values in the x direction and the y direction obtained by loading the field source in the first polarization direction,respectively representing the electric field values in the x direction and the y direction obtained by loading the field source in the second polarization direction,respectively representing the magnetic field values in the x direction and the y direction obtained by loading the field source in the first polarization direction,respectively representing the values of the magnetic field in the x and y directions obtained by loading the field source in the second polarization direction. The expression of the apparent resistivity is as follows,

in the formula (I), the compound is shown in the specification,denotes the apparent resistivity in the ij direction, ω denotes the angular frequency, μ₀＝4π×10^-7H/m；

In the formula (I), the compound is shown in the specification,

indicating the phase response in the ij direction.

Based on the above principle, the magnetotelluric forward modeling method based on the geospatial solution strategy provided by the embodiment of the invention comprises the following steps:

step 1: according to the DEM data, constructing a geoelectric geometric model of the three-dimensional MT in the solution area, and establishing physical parameters of different grid areas;

the physical parameters are resistivity parameters, dielectric constant parameters and permeability parameters of the grid region, and the DEM data are geodetic satellites, InSAR data and topographic maps to obtain Digital Evaluation Map (DEM) data of the measured region.

Step 2: the geoelectric geometric model was discretized into a series of non-overlapping using the Tetgen program. Wherein each tetrahedral unit is attached to a region in the three-dimensional MT geoelectric geometric model.

And step 3: acquiring an integral weak solution form satisfied by A-phi-psi coupling potential three-dimensional MT forward modeling based on a geospatial solution strategy, and performing node-based discrete conversion on the integral weak solution form by adopting a node basis function and a node of a tetrahedral unit to obtain a sparse linear equation set;

and 4, step 4: loading incident fields in two different polarization directions, placing the incident fields in a right-end item of a sparse linear equation set, and solving the sparse linear equation to obtain a magnetic field vector position A, an electric field scalar position phi and a magnetic field scalar position psi which correspond to each node of a tetrahedron unit in a solution area in the two different polarization directions;

and 5: and calculating the combination of one or more parameters of the electric field, the magnetic field, the impedance tensor, the apparent resistivity and the phase response of the observation point according to the magnetic field vector position A, the electric field scalar position phi and the magnetic field scalar position psi.

In addition, the invention also provides a device for realizing the method, which comprises the following steps: the system comprises a model construction module, a tetrahedral unit dispersion module and an operation module, wherein the model construction module is used for constructing a geoelectricity geometric model for solving the three-dimensional MT of the region; the tetrahedral unit discretization module is used for discretizing the geoelectrical geometric model into a series of non-overlapping tetrahedral units; the operation module is used for performing forward calculation based on the A-phi-psi coupling potential of the geospatial solution strategy to finally obtain the target parameters of the observation points. The present invention is not limited to the specific contents of the method.

The invention also provides a computer device which comprises a memory and a processor, wherein the memory is stored with a computer program capable of running on the processor, and the processor executes the computer program to realize the magnetotelluric forward modeling method based on the geospatial solution strategy, the invention also provides a storage medium which is stored with the computer program, and the computer program is executed by the processor to execute the steps of the magnetotelluric forward modeling method based on the geospatial solution strategy, in order to further verify the accuracy and convergence of the developed system for solving the topographic problem, a topographic model like a trapezoidal mountain in figure 2 is established, a trapezoidal topographic model is arranged at the center of coordinate axes, the height of the trapezoid is 450m, the width of the top of the trapezoid is 450m × 450m, the width of the bottom of the trapezoid is 2000m × 2000m, the size of the whole solving area is 20km × 20km × 20km, and the conductivity of the pure topographic model is 0.01 s.m^-1The electrical conductivity of air is 1e^-8s·m^-120 observation sections are arranged along the y axis, and the range of each observation section is y ∈ [ -1500m,1500m]Excitation frequencies of 2Hz, in particular in the x-y and x-z sections of FIG. 3As shown, fig. 4 shows a mesh generation diagram. As can be seen from FIG. 5, apparent resistivity ρ_xxThe numerical value is very small, and the high value area represents the position corresponding to the ridge; and apparent resistivity ρ_xyAnd low resistance phenomenon is shown at the position of the ridge. Apparent resistivity ρ_xyThe interface in the x direction is obvious, and the interface in the y direction is difficult to have one-to-one correspondence with the actual topographic relief in the apparent resistivity contour plane diagram. As can be seen from FIG. 6, apparent resistivity ρ_yyThe numerical value is very small, and the high value area represents the position corresponding to the ridge; apparent resistivity ρ_yxThe contour map is apparent at the y-direction interface.

It should be emphasized that the examples described herein are illustrative and not restrictive, and thus the invention is not to be limited to the examples described herein, but rather to other embodiments that may be devised by those skilled in the art based on the teachings herein, and that various modifications, alterations, and substitutions are possible without departing from the spirit and scope of the present invention.