阅读说明：本技术 一种三维大地电磁正演数值模拟方法 (Three-dimensional magnetotelluric forward modeling numerical simulation method ) 是由 杨刚强 柳建新 郭荣文 王永斐 于 2021-09-22 设计创作，主要内容包括：本申请公开了一种三维大地电磁正演数值模拟方法,首先根据目标地质体构建电阻率分布模型,然后通过多重网格法粗化电阻率分布模型,并通过交错网格有限差分法离散粗细网格上的双旋度方程,得到系数矩阵,接着对离散后的双旋度方程作4色分块处理,通过二维有限差分法计算出目标地质体模型的边界条件,并计算出细网格上双旋度方程的右端项,最后使用基于4色分块高斯-赛德平滑的多重网格法求解；更换极化模式并重复上述过程,根据不同极化模式下的电场分量和磁场分量,计算出对应测点的视电阻率和阻抗相位。本申请采用线分块高斯-赛德平滑技术,有效去除粗细网格上的高频残差分量,提高了多重网格法的计算效率,达到快速收敛的目的。(The application discloses a three-dimensional magnetotelluric forward modeling numerical simulation method, which comprises the steps of firstly constructing a resistivity distribution model according to a target geologic body, then coarsening the resistivity distribution model by a multi-grid method, dispersing a double-rotation equation on a coarse grid and a fine grid by a staggered grid finite difference method to obtain a coefficient matrix, then carrying out 4-color blocking processing on the dispersed double-rotation equation, calculating the boundary condition of the target geologic body model by a two-dimensional finite difference method, calculating the right-end term of the double-rotation equation on the fine grid, and finally solving by using a multi-grid method based on 4-color blocking Gaussian-Saidel smoothing (ii) a Changing the polarization mode and repeating the above process according to the electric field in different polarization modesMeasuring the magnetic field component, and calculating apparent resistivity and impedance phase of the corresponding measuring point. The method adopts the line blocking Gaussian-Saider smoothing technology, effectively removes the high-frequency residual error component on the coarse and fine grids, improves the calculation efficiency of the multiple grid method, and achieves the purpose of rapid convergence.)

1. A three-dimensional magnetotelluric forward modeling numerical simulation method is characterized by comprising the following steps:

step 1) constructing a corresponding resistivity distribution model according to the shape, the size and the resistivity distribution characteristics of a target geologic body;

step 2) coarsening the resistivity distribution model by a multi-grid method to obtain a resistivity distribution model of a coarse grid and a fine grid under different scales, and calculating a limiting operator and an interpolation operator of a spatial relationship between the fine grid and the coarse grid;

step 3) dispersing a double rotation equation satisfied by an electric field on a thick grid and a thin grid by a staggered grid finite difference method, and performing 4-color blocking processing on the dispersed double rotation equation based on a line blocking Gaussian-Saider smoother principle;

step 4) calculating the boundary condition of the target geological body model by a two-dimensional finite difference method according to the resistivity distribution model of the fine grid, and calculating the right-end term of a double rotation equation on the fine grid;

step 5) according to the double rotation equation after the dispersion on the fine grid, performing pre-smoothing by a line blocking Gaussian-Saider smoother to obtain the electric field component on the corresponding fine grid, and calculating the electric field residual component on the fine grid;

step 6) calculating the electric field residual error component on the coarse grid according to the electric field residual error component on the fine grid and the limiting operator;

step 7) solving the double rotation equation according to the double rotation equation after the dispersion on the coarse grid and the electric field residual error component to obtain an electric field iteration quantity on the coarse grid;

step 8) calculating the electric field iteration quantity on the fine grid according to the electric field iteration quantity on the coarse grid and the interpolation operator;

step 9) calculating a new electric field component according to the electric field component and the electric field iteration quantity on the fine grid, smoothing the electric field component by using a line blocking Gaussian-Saider smoother, and updating the electric field component on the fine grid;

step 10) repeating the step 5) to the step 9) until the relative residual error on the fine grid is smaller than a convergence threshold value;

step 11) calculating the magnetic field component on the fine grid according to the electric field component and the rotation operator on the fine grid;

step 12) changing the polarization mode, and repeating the step 4) to the step 11);

and step 13) calculating apparent resistivity and impedance phase of the corresponding measuring point according to the electric field component and the magnetic field component under different polarization modes on the fine grid, thereby completing three-dimensional magnetotelluric forward numerical simulation of a calculated frequency.

2. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 1, wherein the step 1) comprises the following specific processes:

determining a research area, calculating frequency and an observation mode, arranging a corresponding measuring network and constructing a three-dimensional cuboid model;

the three-dimensional cuboid model is divided into a plurality of slender cube units which are respectively arranged along、、The number of elongated cube elements of the directional subdivision is shown as、Andthe length, width and height of each slender square unit are respectively、And；

and assigning the resistivity of each elongated square body unit according to the shape, the size and the resistivity distribution characteristics of the target geologic body, and further realizing the construction of a resistivity distribution model corresponding to the target geologic body.

3. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 2, wherein the step 2) comprises the following specific processes:

will be provided with、、The elongated square units in three directions are combined in pairs to form a thick rectangular unit, namely the eight elongated square units are combined into a thick rectangular unit and then are arranged along、、The number of the thick cuboid units which are divided in the direction is、And；

the length, width, height and resistivity of the thick cuboid units are determined by the length, width, height and resistivity of the thin and long cuboid units contained in the thick cuboid units, and then a resistivity distribution model of the thick grid is constructed;

respectively calculating limit operators according to the spatial relationship between the elongated square units and the thick rectangular unitsAnd interpolation operator；

Based on the requirement of three-dimensional magnetotelluric forward numerical simulation, the model can be further coarsened to obtain a coarser resistivity distribution model and obtain a corresponding limiting operator、……And interpolation operator、……，The number of model coarsening times.

4. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 3, wherein the discrete double rotation equation in the step 3) comprises the following specific processes:

the double rotation equation satisfied by the three-dimensional magnetotelluric electric field intensity can be derived from the Max equation:

(1)

in the formula:for the strength of the electric field,is an angular frequency and is associated with the calculated frequencyIn a relationship of，In order to have a magnetic permeability,is the conductivity and resistivity of the mediumIn a relationship of；

With the staggered mesh finite difference method, equation (1) can be discretized into:

(2)

(3)

in the formula:is the electric field component of the middle point of the upper edge of the cuboid unit,is a rotation operator from the middle point of the edge to the center of the surface,is the rotation operator from the center of the face to the midpoint of the edge,the conductivity volume of the cuboid monomer with connected edges is average;

representing electric field components on six faces of the three-dimensional rectangular solid model constructed in step 1) asThat is, the electric field component located inside the three-dimensional rectangular solid model, i.e., the boundary electric field componentExpressed as an internal electric field component, equation (2) can be expressed as:

(4)

in the formula:andare respectively asFrom the coefficient matrix after separation of the internal and boundary electric field componentsItem replacement byAnd moves to the right of the equation as the right term of the equation, equation (4) reduces to:

(5)

the double rotation degree equation satisfied by the electric field on the coarse and fine grids is dispersed by the method to obtain the corresponding coefficient matrix, and the coefficient matrix is used、Andrespectively representing coefficient matrix, electric field component and right-hand term on fine grid、……Representing a matrix of coefficients on a coarse grid.

5. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 4, wherein the 4-color blocking processing on the discretized bi-rotation equation in the step 3) comprises the following specific processes:

based on the working principle of the line blocking Gauss-Saider smoother, the obtained electric field components are divided into 4 types and are marked by different colors, the electric field components connected with all nodes on the same color line are mutually independent and are synchronously calculated, and the boundary coupling is carried out between different colors through the electric field components connected with all nodes on the adjacent lines, namely

Will be provided withDecomposed into 4 small coefficient matrices、、Andwill beIs decomposed into、、Andwill beIs decomposed into、、Andthe solution for equation (5) goes to the solution of the following 4 sub-equations:

(6)。

6. the three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 5, wherein the step 4) comprises the following specific procedures:

selectingOrPolarizing, namely dividing the three-dimensional cuboid model into a series of two-dimensional cuboid models;

calculating the electric field components of the two-dimensional cuboid model by a two-dimensional finite difference method, further obtaining all the electric field components after the three-dimensional cuboid model is divided, and assigning the electric field components positioned on six surfaces of the three-dimensional cuboid model to the electric field componentsAccording toCalculating the right-hand term of the double rotation degree equation, and assigning the electric field component in the inner part to the right-hand termAs the initial estimation valueAnd (4) showing.

7. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 6, wherein the steps 5) -9) comprise the following specific processes:

will be provided withSubstituting into formula (5), and smoothly updating with formula (6) to obtain electric field component on the corresponding fine gridAnd simultaneously calculating the electric field residual components on the fine grid:

(7)

coordination of constraint operatorsThe electric field residual components on the fine grid are limited to those on the coarse grid:

(8)

coefficient matrix obtained from the dispersion on the coarse gridCalculating the electric field iteration quantity on the coarse grid by using a direct solver:

(9)

fitting interpolation operatorsAnd extending the electric field iteration quantity on the coarse grid into the electric field iteration quantity on the fine grid:

(10)

according to electric field component on fine gridAnd amount of electric field iterationCarrying out coarse grid correction on the electric field component to obtain a new electric field component;

substituting the corrected new electric field component into formula (5), and smoothly updating the new electric field component with formula (6) to obtain new electric field component on the corresponding fine grid。

8. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 7,in said step 10), when the electric field residual components on the fine mesh satisfy the relative residual＜1×10^-10And stopping iterative computation.

9. The three-dimensional magnetotelluric forward modeling numerical simulation method according to claim 8, wherein the step 13) comprises the following specific procedures:

according to the electric field component and the magnetic field component under two polarization modes, calculating、Two modes of impedance:

(11)

in the formula:、、andis composed ofCalculation of polarization、The directional electric field component and the magnetic field component,、、andis composed ofCalculation of polarization、Directional electric and magnetic field components;

calculated from the obtained impedance、Apparent resistivity and impedance phase in two modes:

(12)。

Technical Field

The invention relates to the technical field of geophysical technology, in particular to a three-dimensional magnetotelluric forward modeling numerical simulation method.

Background

The geoelectromagnetic prospecting is a method of exploiting the widely distributed naturally varying electromagnetic field (10) in the earth^-4~10⁴Hz), has the advantages of large detection depth, high working efficiency, low expenditure cost and the like, is widely applied to the fields of engineering exploration, metal and oil gas resource exploration, deep geodynamics research and the like, and can efficiently carry out electromagnetic three-dimensional forward inversion, which is a premise that the method can be widely applied.

In the process of forward numerical simulation of the three-dimensional electromagnetic method, a double-rotation equation related to the magnetic field rotation and the electric field rotation needs to be solved, but the equation has abundant null space, and when the frequency is very low, the equation cannot simulate the electrical change, so that the conventional iteration method such as BICGSab is adopted to solve the equation, the convergence efficiency is low, and divergence correction needs to be carried out to accelerate the convergence speed. However, as the calculation frequency is lower and the number of the subdivision grids is more, the iteration number shows linear and rapid increase, so that the conventional iteration method is used for simulating and researching the ultra-deep and complex geological abnormal body, and the high-efficiency production efficiency is difficult to achieve.

The prior document discloses the application of a multiple-grid method in a diffusion electromagnetic field, lays a foundation for the application of the multiple-grid method in magnetotelluric simulation, and considers that the convergence rate of the multiple-grid method does not show linear increase along with the increase of the number of grid nodes, so that the multiple-grid method is more suitable for solving a large equation set compared with a conventional iteration method. But the computational efficiency of the method depends on a smoothing algorithm to a great extent, and the computational efficiency of the conventional smoothing algorithm is low for the double rotation equation of the electromagnetic field.

In summary, researchers need to perform targeted optimization on the existing simulation method, and solve the problems that the calculation scale is increased, the calculation efficiency is reduced in a nonlinear manner, the calculation frequency is lower, the problem is solved more ill, the convergence speed of the conventional iteration method is slow, and the like, so that the requirement of rapid inversion imaging on complex geologic bodies in the exploration of the magnetotelluric method is met.

Disclosure of Invention

The application aims to provide a three-dimensional magnetotelluric forward modeling numerical simulation method, which can ensure result precision and calculation efficiency when electric field and magnetic field response generated by a deep complex geologic body is simulated. The technical scheme of the application is as follows:

a three-dimensional magnetotelluric forward modeling numerical simulation method comprises the following steps:

step 1) constructing a corresponding resistivity distribution model according to the shape, the size and the resistivity distribution characteristics of a target geologic body;

step 2) coarsening the resistivity distribution model by a multi-grid method to obtain a resistivity distribution model of a coarse grid and a fine grid under different scales, and calculating a limiting operator and an interpolation operator of a spatial relationship between the fine grid and the coarse grid;

step 3) dispersing a double rotation equation satisfied by an electric field on a thick grid and a thin grid by a staggered grid finite difference method, and performing 4-color blocking processing on the dispersed double rotation equation based on a line blocking Gaussian-Saider smoother principle;

step 4) calculating the boundary condition of the target geological body model by a two-dimensional finite difference method according to the resistivity distribution model of the fine grid, and calculating the right-end term of a double rotation equation on the fine grid;

step 5) according to the double rotation equation after the dispersion on the fine grid, performing pre-smoothing by a line blocking Gaussian-Saider smoother to obtain the electric field component on the corresponding fine grid, and calculating the electric field residual component on the fine grid;

step 6) calculating the electric field residual error component on the coarse grid according to the electric field residual error component on the fine grid and the limiting operator;

step 7) solving the double rotation equation according to the double rotation equation after the dispersion on the coarse grid and the electric field residual error component to obtain an electric field iteration quantity on the coarse grid;

step 8) calculating the electric field iteration quantity on the fine grid according to the electric field iteration quantity on the coarse grid and the interpolation operator;

step 9) calculating a new electric field component according to the electric field component and the electric field iteration quantity on the fine grid, smoothing the electric field component by using a line blocking Gaussian-Saider smoother, and updating the electric field component on the fine grid;

step 10) repeating the step 5) to the step 9) until the relative residual error on the fine grid is smaller than a convergence threshold value;

step 11) calculating the magnetic field component on the fine grid according to the electric field component and the rotation operator on the fine grid;

step 12) changing the polarization mode, and repeating the step 4) to the step 11);

and step 13) calculating apparent resistivity and impedance phase of the corresponding measuring point according to the electric field component and the magnetic field component under different polarization modes on the fine grid, thereby completing three-dimensional magnetotelluric forward numerical simulation of a calculated frequency.

In some specific embodiments, the step 1) includes the following specific processes:

determining a research area, calculating frequency and an observation mode, arranging a corresponding measuring network and constructing a three-dimensional cuboid model;

the three-dimensional cuboid model is divided into a plurality of slender cube units which are respectively arranged along、、The number of elongated cube elements of the directional subdivision is shown as、Andthe length, width and height of each slender square unit are respectively、And；

and assigning the resistivity of each elongated square body unit according to the shape, the size and the resistivity distribution characteristics of the target geologic body, and further realizing the construction of a resistivity distribution model corresponding to the target geologic body.

In some specific embodiments, the step 2) includes the following specific processes:

will be provided with、、The elongated square units in three directions are combined in pairs to form a thick rectangular unit, namely the eight elongated square units are combined into a thick rectangular unit and then are arranged along、、The number of the thick cuboid units which are divided in the direction is、And；

the length, width, height and resistivity of the thick cuboid units are determined by the length, width, height and resistivity of the thin and long cuboid units contained in the thick cuboid units, and then a resistivity distribution model of the thick grid is constructed;

respectively calculating limit operators according to the spatial relationship between the elongated square units and the thick rectangular unitsAnd interpolation operator；

Based on the requirement of three-dimensional magnetotelluric forward numerical simulation, the model can be further coarsened to obtain a coarser resistivity distribution model and obtain a corresponding limiting operator、……And interpolation operator、……，The number of model coarsening times.

In some specific embodiments, the discrete bi-rotation equation in step 3) includes the following specific processes:

the double rotation equation satisfied by the three-dimensional magnetotelluric electric field intensity can be derived from the Max equation:

(1)

in the formula:for the strength of the electric field,is an angular frequency and is associated with the calculated frequencyIn a relationship of，In order to have a magnetic permeability,is the conductivity and resistivity of the mediumIn a relationship of；

With the staggered mesh finite difference method, equation (1) can be discretized into:

(2)

(3)

in the formula:is the electric field component of the middle point of the upper edge of the cuboid unit,is a rotation operator from the middle point of the edge to the center of the surface,is the rotation operator from the center of the face to the midpoint of the edge,the conductivity volume of the cuboid monomer with connected edges is average;

representing electric field components on six faces of the three-dimensional rectangular solid model constructed in step 1) asThat is, the electric field component located inside the three-dimensional rectangular solid model, i.e., the boundary electric field componentExpressed as an internal electric field component, equation (2) can be expressed as:

(4)

in the formula:andare respectively asFrom the coefficient matrix after separation of the internal and boundary electric field componentsItem replacement byAnd moves to the right of the equation as the right term of the equation, equation (4) reduces to:

(5)

the double rotation equation satisfied by electric field on coarse and fine grids is dispersed by the method to obtain corresponding coefficient matrix, and the coefficient matrix is used、Andrespectively representing coefficient matrix, electric field component and right-hand term on fine grid、……Representing a matrix of coefficients on a coarse grid.

In some specific embodiments, the 4-color blocking processing on the discretized bi-rotation equation in step 3) includes the following specific processes:

based on the working principle of the line blocking Gauss-Saider smoother, the obtained electric field components are divided into 4 types and are marked by different colors, the electric field components connected with all nodes on the same color line are mutually independent and are synchronously calculated, and the boundary coupling is carried out between different colors through the electric field components connected with all nodes on the adjacent lines, namely

Will be provided withDecomposed into 4 small coefficient matrices、、Andwill beIs decomposed into、、Andwill beIs decomposed into、、Andthe solution for equation (5) goes to the solution of the following 4 sub-equations:

(6)。

in some specific embodiments, the step 4) includes the following specific processes:

selectingOrPolarizing, namely dividing the three-dimensional cuboid model into a series of two-dimensional cuboid models;

calculating the electric field components of the two-dimensional cuboid model by a two-dimensional finite difference method, further obtaining all the electric field components after the three-dimensional cuboid model is divided, and assigning the electric field components positioned on six surfaces of the three-dimensional cuboid model to the electric field componentsAccording toCalculating the right-hand term of the double rotation degree equation, and assigning the electric field component in the inner part to the right-hand termAs the initial estimation valueAnd (4) showing.

In some specific embodiments, the steps 5) to 9) include the following specific processes:

will be provided withSubstituting into formula (5), and smoothly updating with formula (6) to obtain electric field component on the corresponding fine gridAnd simultaneously calculating the electric field residual components on the fine grid:

(7)

coordination of constraint operatorsThe electric field residual components on the fine grid are limited to those on the coarse grid:

(8)

coefficient matrix obtained from the dispersion on the coarse gridCalculating the electric field iteration quantity on the coarse grid by using a direct solver:

(9)

fitting interpolation operatorsAnd extending the electric field iteration quantity on the coarse grid into the electric field iteration quantity on the fine grid:

(10)

according to electric field component on fine gridAnd amount of electric field iterationCarrying out coarse grid correction on the electric field component to obtain a new electric field component;

substituting the corrected new electric field component into formula (5), and smoothly updating the new electric field component with formula (6) to obtain new electric field component on the corresponding fine grid。

In some specific embodiments, in the step 10), when the electric field residual components on the fine mesh satisfy the relative residual＜1×10^-10And stopping iterative computation.

In some specific embodiments, the step 13) includes the following specific processes:

according to the electric field component and the magnetic field component under two polarization modes, calculating、Two modes of impedance:

(11)

in the formula:、、andis composed ofCalculation of polarization、The directional electric field component and the magnetic field component,、、andis composed ofCalculation of polarization、Directional electric and magnetic field components;

calculated from the obtained impedance、Apparent resistivity and impedance phase in two modes:

(12)。

the technical scheme provided by the application has at least the following beneficial effects:

1. the method adopts a line blocking Gaussian-Saider smoothing technology with high parallelism, so that the method contains information that the divergence of electric field components in magnetotelluric digital simulation is zero, divergence correction is not needed, and the calculation efficiency of multiple grids is accelerated; in addition, high-frequency residual components on the coarse and fine grids can be effectively removed, the coefficient matrixes of the double-rotation equation are subjected to reordering decomposition by adopting a line blocking Gauss-Saider smoother, and the sparsity of the decomposed 4 small coefficient matrixes is improved, so that the method has higher convergence speed and calculation efficiency at low frequency.

2. According to the method, the cuboid units are used for subdividing the three-dimensional geologic body model, and the resistivity assigned to each cuboid unit can be different, so that the target geologic body with any resistivity distribution can be simulated.

3. In the present application: the double rotation degree equation satisfied by the electric field on the cuboid unit is dispersed by using the staggered grid finite difference method, which is beneficial to simply and efficiently dispersing the double rotation degree equation on the thick and thin grid, and further obtaining the corresponding coefficient matrix; a two-dimensional finite difference numerical solution is used as a boundary condition, so that a numerical simulation result is more accurate; the dual-rotation equation is obtained by solving the finite difference method dispersion of the staggered grids by using the multiple-grid method, so that the low-frequency electric field residual components on the fine grid can be limited to the high-frequency residual components on the coarse grid for elimination, and the purpose of rapid convergence is achieved.

Drawings

In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings that are needed in the description of the embodiments or the prior art will be briefly described below, it being understood that the following drawings only illustrate some embodiments of the present application and therefore should not be considered as limiting the scope, and that other drawings may be derived from those drawings without inventive effort for a person skilled in the art.

Fig. 1 is a flowchart of a three-dimensional magnetotelluric forward numerical simulation method provided in the present application;

FIG. 2 is a block diagram of a line-blocking Gaussian-Sade smoother in the method of the present application, wherein FIG. 2 (a) is a line smoothing system, FIG. 2 (b) is a schematic diagram of line blocking, and 1, 2, 3, and 4 represent 4 different color marks, respectively;

FIG. 3 is a schematic diagram of a low-high resistance combined anomaly model in example 1 of the present application;

FIG. 4 is the result of the low-high impedance combined anomaly on the surface survey line calculated by the rapid multigrid method and the International publication program ModEM in example 1 of the present application, wherein FIG. 4 (a) isThe comparison of the apparent resistivities obtained in the two modes, FIG. 4 (b)The phase of the two impedance phases obtained in the mode are compared with each other, and FIG. 4 (c) is a graphThe comparison of the apparent resistivities obtained in the two modes is shown in FIG. 4 (d)Comparing the two impedance phases in the mode;

fig. 5 is a graph comparing the convergence curves of the BICGstab with and without divergence correction in example 1 of the present application.

Detailed Description

In order to facilitate understanding of the present application, the technical solutions in the present application will be described more fully and in detail with reference to the drawings and the preferred embodiments, but the scope of protection of the present application is not limited to the following specific embodiments, and all other embodiments obtained by a person of ordinary skill in the art based on the embodiments in the present application without creative efforts shall fall within the scope of protection of the present application.

It will be understood that when an element is referred to as being "coupled" or "connected" to another element, it can be directly coupled, connected or communicated with the other element or indirectly coupled, connected or communicated with the other element via other intervening elements.

Unless otherwise defined, all terms of art used hereinafter have the same meaning as commonly understood by one of ordinary skill in the art. The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the scope of the present application.

Referring to fig. 1, a three-dimensional magnetotelluric forward modeling numerical simulation method specifically comprises the following steps:

firstly, constructing a three-dimensional resistivity distribution model

Determining a research area, calculating frequency and an observation mode according to the shape, the size and the resistivity distribution of the target geologic body, arranging a corresponding measuring net and constructing a three-dimensional cuboid model.

The three-dimensional rectangular solid model is subjected to mesh subdivision to form a plurality of slender square units which are respectively arranged along the edge、、The number of elongated cube elements of the directional subdivision is shown as、Andthe length, width and height of each slender square unit are respectively、And。

the resistivity of each elongated cube unit is assigned according to the shape, the size and the resistivity distribution characteristics of the target geologic body, the resistivity of each elongated cube unit can be different, so that the target geologic body under any condition can be simulated to construct a corresponding resistivity distribution model, and the resistivity of each elongated cube unit in the air layer is assigned to be。

Second step, coarsening three-dimensional resistivity distribution model

Will be provided with、、The elongated square units in three directions are combined in pairs to form a thick rectangular unit, namely the eight elongated square units are combined into a thick rectangular unit and then are arranged along、、The number of the thick cuboid units which are divided in the direction is、And。

the length, width and height of the thick cuboid units are obtained by correspondingly adding the length, width and height of the contained slender cuboid units, and the resistivity of the thick cuboid units is obtained by averaging the volume of the resistivity of the contained slender cuboid units, so that a resistivity distribution model of the thick grid is constructed.

Respectively calculating limit operators according to the spatial relationship between the elongated square units and the thick rectangular unitsAnd interpolation operator。

Based on the requirement of three-dimensional magnetotelluric forward numerical simulation, the grid can be coarsened further to obtain a coarser resistivity distribution model and obtain a corresponding limiting operator、……And interpolation operator、……，The number of model coarsening times.

Thirdly, a bispin equation on a discrete coarse-fine grid

The double rotation equation satisfied by the three-dimensional magnetotelluric electric field intensity can be derived from the Max equation:

(1)

in the formula:for the electric field strength (V/m),is an angular frequency (Hz) and is related to the calculated frequencyIn a relationship of，Magnetic permeability (H/m), without considering underground rock magnetization，Is dielectric conductivity (S/m) and is related to resistivity（) In a relationship of。

With the staggered mesh finite difference method, equation (1) can be discretized into:

(2)

(3)

in the formula:is the electric field component of the middle point of the upper edge of the cuboid unit,is a rotation operator from the middle point of the edge to the center of the surface,is the rotation operator from the center of the face to the midpoint of the edge,the conductivity volume of the cuboid monomer connected with the edge is average.

Representing electric field components on six faces of the constructed three-dimensional rectangular solid model asThat is, the electric field component located inside the three-dimensional rectangular solid model, i.e., the boundary electric field componentExpressed as an internal electric field component, equation (2) can be expressed as:

(4)

in the formula:andare respectively asFrom the coefficient matrix after separation of the internal and boundary electric field componentsItem replacement byAnd moves to the right of the equation as the right term of the equation, equation (4) reduces to:

(5)

the double rotation degree equation satisfied by the electric field on the coarse and fine grids is dispersed by the method to obtain the corresponding coefficient matrix, and the coefficient matrix is used、Andrespectively representing coefficient matrix, electric field component and right-hand term on fine grid、……Representing a matrix of coefficients on a coarse grid.

Fourth step, line block Gaussian-Saider smoother operating principle deduction

The strength of the electric field in the magnetotelluric can be deduced from the Maxwell equation setAnd the intensity of the magnetic fieldThe relationship of (1):

(13)

as can be seen from equation (13), the electric field strength can be obtained from the curl of the magnetic field strength, and the magnetic field strength can be obtained from the curl of the electric field strength.

The scattered electric field component is subjected to loop integration along the surface of a cuboid unit to obtain a magnetic field component of the element surface central point; and the magnetic field component is subjected to loop integration along double volume units in the staggered grid finite difference method to obtain the electric field component of the middle point of the edge of the cuboid unit.

Thus, is obtained byThe electric field component connecting all nodes on one line in the direction is only related to the electric field component of the nodes on the adjacent line, as shown in FIG. 2a, the systemThe arrow mark on the periphery represents the electric field component on the adjacent line, and the arrow mark inside the system and connected with the center represents the electric field component connected with all nodes on the current line to be solved.

The electric field components connected to all nodes on two lines that are not adjacent to each other are independent of each other, and therefore, the electric field components in the third step can be classified into 4 types, as indicated by 1, 2, 3, and 4 in fig. 2b, and the same symbols are represented as the same type and shown in one color. The electric field components connected with all nodes on the same color line are mutually independent and can be synchronously calculated, and the electric field components connected with all nodes on the adjacent lines among different colors are subjected to boundary coupling.

Then, can beDecomposed into 4 small coefficient matrices、、Andwill beIs decomposed into、、Andwill beIs decomposed into、、Andcolor mapping of No. 1Color mapping of No. 2Color mapping of No. 3Color mapping of No. 4. The solution for equation (5) is shifted to a solution of the following 4 sub-equations:

(6)

the boundary coupling is performed on the electric field components connected by all nodes on the adjacent lines among the different colors, namely, the equation in the formula (6) is sequentially calculated: first from the initial valueIs decomposed outCalculated according to the first equation of equation (6)Update(ii) a Then from the updatedIs decomposed to obtainCalculated according to the second expression in equation (6)Update(ii) a Then from after updatingIs decomposed to obtainCalculated according to the third equation in equation (6)Update(ii) a Finally from updatedIs decomposed to obtainCalculated according to the last expression in equation (6)UpdateI.e. a line blocking gaussian-seidel smoothing is completed.

In order to prevent the repeated calculation of LU decomposition of coefficient matrix in the calculation process, the method is to calculate in advance、、Andthe LU decomposition result of (1) is saved.

Fifthly, calculating the right-end term of the discrete double rotation equation

Cutting the three-dimensional rectangular model into a series of two-dimensional rectangular models according to the polarization mode, wherein,the polarization is alongCutting direction layer by layerA two-dimensional rectangular model of a plane,the polarization is alongCutting direction layer by layerPlanar two-dimensional rectangleAnd (4) modeling.

Calculating the electric field components of the two-dimensional cuboid model by a two-dimensional finite difference method, further obtaining all the electric field components after the three-dimensional cuboid model is divided, and assigning the electric field components positioned on six surfaces of the three-dimensional cuboid model to the electric field componentsAccording toCalculating the right-hand term of the double rotation degree equation, and assigning the electric field component in the inner part to the right-hand termAs the initial estimation valueAnd (4) showing.

Sixth step, solving formula (5)

The multi-grid method limits the low-frequency residual components on the fine grid to the high-frequency residual components on the coarse grid for removing, and the line blocking Gaussian-Sade smoother can quickly and effectively remove the high-frequency residual components on the grid. Equation (5) is therefore solved using a fast multiple-grid method based on a line-blocking gaussian-seidel smoother.

Pre-smoothing using a line-blocking Gaussian-Sade smoother, i.e. in the fifth stepSubstituting into formula (5), and smoothly updating with formula (6) to obtain electric field component on the corresponding fine gridAnd simultaneously calculating the electric field residual components on the fine grid:

(7)

coordination of constraint operatorsThe electric field residual components on the fine grid are limited to those on the coarse grid:

(8)

coefficient matrix obtained from the dispersion on the coarse gridCalculating the electric field iteration quantity on the coarse grid by using a direct solver:

(9)

fitting interpolation operatorsAnd extending the electric field iteration quantity on the coarse grid into the electric field iteration quantity on the fine grid:

(10)

according to electric field component on fine gridAnd amount of electric field iterationCoarse grid correction of electric field componentsObtaining new electric field components;

using a line blocking Gaussian-Saider smoother to advanceAnd (4) smoothing after line, namely substituting the corrected new electric field component into a formula (5), and smoothly updating the new electric field component by using a formula (6) to obtain a new electric field component on the corresponding fine grid。

Step seven, iterative cycle

Repeating the calculation process of the sixth step until the electric field residual error components on the fine mesh satisfy the relative residual errors＜1×10^-10Stopping the iterative calculation, and calculating the electric field component according to the formula (13)The rotation degree is taken to obtain the magnetic field component in the geologic body model。

In order to prevent the direct solver from repeatedly calculating the LU decomposition of the coefficient matrix, the LU decomposition result of the coefficient matrix is stored in the computer in advance.

And eighthly, replacing the polarization mode, and repeating the fifth step, the sixth step and the seventh step.

Ninth, calculating apparent resistivity and impedance phase of measuring point on ground corresponding measuring line

According to the electric field component and the magnetic field component under two polarization modes, calculating、Two modes of impedance:

(11)

in the formula:、、andis composed ofCalculation of polarization、The directional electric field component and the magnetic field component,、、andis composed ofCalculation of polarization、Directional electric field components and magnetic field components.

Calculated from the obtained impedance、Apparent resistivity and impedance phase in two modes:

(12)

the three-dimensional magnetotelluric forward modeling numerical simulation method provided by the application is examined below.

Example 1

In order to examine the method proposed by the present invention, a low-high resistance combined anomaly model as shown in fig. 3 was designed, the details of which are as follows: with low-resistance abnormal body and high-resistance abnormal body edge in the groundThe sizes of the combined abnormal bodies which are transversely arranged in the direction are 10km multiplied by 10km, the buried depths of the combined abnormal bodies are 10km, and the distance between the combined abnormal bodies is 10 km; the air layer resistivity of the three-dimensional low-high resistance combined abnormal body model isBackground resistivity of 100 in large earth formationThe resistivity of the low-resistance abnormal body is 10The resistivity of the high-resistance abnormal body is 1000. Taking the projection point of the midpoint between the two abnormal bodies on the ground as the seatThe standard origin point O and the measuring line are，And 25 measuring points are uniformly arranged at intervals within the range of-25.5 km to 25.5km in direction.

If the ground layer is subdivided by using 32 × 32 × 32 rectangular parallelepiped units with the size of 2.5km × 2.5km × 2.5km, the simulation area is:andthe directions are all-40 km to 40km,the direction is 0km to 80 km; if the ground layer is divided by using 64 × 64 × 64 rectangular parallelepiped units with the size of 10km × 10km × 10km, the simulation area is:andthe directions are all-32 km to 32km,the direction is 0km to 64 km; if the ground layer is divided by using 128 × 128 × 128 rectangular parallelepiped units with the size of 10km × 10km × 10km, the simulation area is as follows:andthe directions are all-64 km to 64km,the direction is 0km to 128 km. The three mesh divisions can effectively simulate the low-high resistance combined abnormal body model shown in the figure 3, and similar cuboid units are adopted for the air layer.

FIG. 4 shows the three-dimensional magnetotelluric forward modeling numerical simulation method and the low-high resistance combined abnormal body calculated by the international published ModEM program on the earth surface survey line、Apparent resistivity and impedance phase contrast curves in two modes, wherein the size of a subdivision grid is 64 multiplied by 64, and the calculation frequency is 0.1 Hz. As can be seen from FIG. 4, the method of the present application substantially matches the results of the numerical simulation using the International published ModEM program, and the reliability of the method of the present application in the accuracy of the results is verified.

FIG. 5 is a graph comparing the convergence curves of the present application method and the conventional iterative method BICGSab with and without divergence correction, wherein the size of the subdivision grid is 64 × 64 × 64, the calculation frequency is 0.1Hz, and the calculation polarization mode isFrom fig. 5, it can be known that the method of the present application achieves fast convergence within several multiple grid iterations.

TABLE 1

Table 1 shows the statistical data of the iteration times and the computation time of the fast multiple grid method based on the line blocking gaussian-seider smoother in this embodiment under different computation frequencies, where the size of the subdivision grid is 64 × 64 × 64, and the computation polarization mode isFrom table 1, it can be seen that the calculation efficiency of the method of the present application does not become lower as the calculation frequency becomes lowerAnd the calculation speed of the method is about 6 times of that of the traditional iterative method for correcting the divergence in the range of 0.1-0.001 Hz.

TABLE 2

Table 2 is a statistical table of the iteration times and the calculation time of the fast multiple-grid method based on the line-blocking gaussian-seider smoother in this embodiment under different mesh sizes, where the calculation frequency is 0.1Hz, and the calculation polarization mode is. As can be seen from the table 2, the iteration times of the method cannot be increased along with the increase of the size of the subdivision grids, and the calculation advantages are more obvious when the problem of more subdivision grids is simulated.

Therefore, the method has important significance for developing large-scale three-dimensional magnetotelluric forward numerical simulation research of the complex abnormal body, can greatly improve the calculation efficiency, and reduces the cost and expense brought by time.

The above description is only a few examples of the present application and does not limit the scope of the claims of the present application, and it will be apparent to those skilled in the art that various modifications and variations can be made in the present application. Any improvement or equivalent replacement directly or indirectly applicable to other related technical fields within the spirit and principle of the present application by using the contents of the specification and the drawings of the present application shall be included in the protection scope of the present application.