Three-dimensional finite element transient electromagnetic forward modeling calculation method based on iterative solution method

Abstract

The invention relates to the technical field of geophysical exploration, and particularly provides a three-dimensional finite element transient electromagnetic forward modeling calculation method based on an iterative solution method.The method comprises the steps that tetrahedron space discretization is carried out based on a vector finite element method, and an electric field is given to six edges of a space tetrahedron; constructing a back-pushing Euler finite element equation set, and assembling a linear system formed by a large sparse matrix based on the equation set; performing iterative solution on the linear system based on an FGMRES algorithm to obtain an electric field value corresponding to a spatial tetrahedron unit edge; optimizing a spatial tetrahedral mesh based on a quadratic nesting method, and optimizing a time step length based on a third-order backward Euler difference format; a grid pre-encryption method is adopted, the calculation error of iterative solution is reduced, a region decomposition method is used for achieving parallel calculation, and the calculation efficiency is greatly improved. According to the method, the requirements of high precision and high efficiency in three-dimensional transient electromagnetic forward modeling calculation can be met, and an advanced calculation tool is provided for resource exploration in a complex geological environment.

Description

Technical Field

[0001] The present invention relates to the technical field of physical exploration, and in particular, to a three-dimensional finite element transient electromagnetic forward modeling calculation method based on an iterative solution method. Background Art

[0002] The transient electromagnetic method is a commonly used geophysical exploration method. It has relatively low requirements for transmitting and receiving devices and is widely applied in many fields, such as mineral exploration, hydro-environmental investigation, tunnel advanced detection, etc. With the continuous progress of detection instruments and the increasing improvement of geological data, there is a requirement for rapid and stable forward simulation of a more refined large-scale underground environment.

[0003] Currently, there are three common types of three-dimensional transient electromagnetic forward modeling methods. One is to calculate the frequency-domain response of the transient source and then convert it to the time domain. Another is the time-domain algorithm that directly performs iterative calculations through discrete time, and the solution method based on the Krylov rational subspace.

[0004] The frequency-domain algorithm calculates the frequency-domain response and then uses frequency-time conversion methods such as Hankel transform to convert it back to the time domain. It is widely used in the calculation of one-dimensional TEM responses. However, for three-dimensional responses, its requirements for the calculation grid and the number of frequency points are relatively high, and the calculation time is also relatively longer. The time-domain algorithm directly performs iterative solutions of Maxwell's equations by discretizing space and time. According to different time discretization methods, it is divided into explicit calculation and implicit calculation. Explicit time-domain algorithms mostly use the finite difference method for calculation, but they need to follow the CFL stability condition, which greatly limits their application. Currently, the more popular algorithm is to use the implicit time-step algorithm. Due to its unconditional convergence property, it is more free to select the time step in the applications of finite difference, finite volume, and finite element algorithms. In recent years, in addition to the above time-frequency domain solution methods, more research has been done on the TEM forward modeling based on the Krylov subspace method. Compared with the above algorithms, the Krylov subspace method avoids time iteration and directly solves the TEM response of the corresponding time trace by constructing the corresponding Krylov subspace with the time-related matrix exponential function.

[0005] Currently, the forward modeling research problems of three-dimensional TEM mainly focus on more accurate field source modeling and more complex underground environment simulation, and mainly face the following problems:

[0006] I. Response calculation of large-scale geological models: Currently, the commonly used frequency-domain and time-domain methods both require solving large sparse equations in the form of Ax = b. Directly solving the LU decomposition requires a large amount of memory and relatively long calculation time. Relatively speaking, the Krylov strategy based on parallel strategies, iterative solution strategies, and model reduction is more suitable for the response calculation of large models. Generally speaking, there are relatively few calculation methods that can be applied to large-scale model calculations currently.

[0007] II. Full-time response calculation: Currently, the conventional calculation method is through convolution calculation with the current waveform or direct calculation in the time domain. The convolution calculation method is the only way to calculate the full-time response by frequency-domain algorithms. It can achieve good results for simple waveforms such as square waves and sine waves, but it is difficult to accurately calculate complex waveforms. The time-domain calculation method is more direct, but it will increase the required calculation step size and reduce the calculation efficiency.

[0008] III. Response calculation considering induced polarization: Since there is an analytical formula for the dispersion model in the frequency domain, the complex resistivity corresponding to the frequency can be directly replaced in the frequency domain, and the calculation is relatively simple, but the requirements for calculation frequency points and grids are also higher, and it is more time-consuming. The commonly used dispersion models are usually expressed as fractional derivatives in the time domain in the time domain. Currently, there are various approximation methods. In addition, there is also a stretched exponential model established directly from the time domain. Although there are already many algorithms, how to calculate the response of polarized media efficiently and accurately is still a problem to be solved. Summary of the Invention

[0009] In view of this, the present invention proposes a three-dimensional finite element transient electromagnetic forward calculation method based on an iterative solution method to solve the problems existing in the above-mentioned prior art.

[0010] To achieve the above object, the present invention proposes a three-dimensional finite element transient electromagnetic forward calculation method based on an iterative solution method, which is characterized by including:

[0011] Based on the vector finite element method, tetrahedral space discretization is performed, and the electric field is assigned to the six edges of the spatial tetrahedron;

[0012] The backward Euler method is used to discretize the time step, a backward Euler finite element equation set is constructed, and a linear system composed of a large sparse matrix is assembled based on the backward Euler finite element equation set;

[0013] Based on the FGMRES algorithm, the linear system is iteratively solved to obtain the electric field values corresponding to the edges of the spatial tetrahedron elements;

[0014] The spatial tetrahedral mesh is optimized based on the quadratic nesting method, and the time step is optimized based on the third-order backward Euler difference format;

[0015] The grid pre - encryption method is adopted to reduce the computational error of iterative solution, and the domain decomposition method is used to achieve parallel computing.

[0016] Furthermore, the finite - element equation set is as follows:

[0017]

[0018] In the formula, M is the mass matrix, S is the stiffness matrix, J s is the field source, i, j represent matrix or vector elements, t represents time, E j represents the electric - field value on the edge, and n represents the time step after time discretization.

[0019] Furthermore, the linear - system matrix is as follows:

[0020] Ax = b

[0021] In the formula, A represents the linear - system matrix, A = 3M+2ΔtS, x = e n , b = M(4e n-1 -e n-2 ), M is the mass matrix, S is the stiffness matrix, e n represents the electric field after time discretization.

[0022] Furthermore, the process of iteratively solving the linear - system matrix based on the FGMRES algorithm includes:

[0023] Constructing a pre - condition matrix for the linear - system matrix based on the auxiliary - space Maxwell pre - conditioner to reduce the matrix condition number, and converting the linear - system matrix into an equivalent equation set;

[0024] Defining an m - dimensional subspace for the equivalent equation set, and representing any vector at an arbitrary position within the spatial tetrahedron element based on a set of bases in the m - dimensional subspace;

[0025] Constructing an iterative formula based on the linear - representation coefficients, expanding the m - dimensional subspace, constructing subspace bases, and solving the linear - representation coefficients based on the subspace bases;

[0026] Calculating the iterative formula based on the linear - representation coefficients to obtain the vector representation at any position within the spatial tetrahedron element.

[0027] Furthermore, the m - dimensional subspace is as follows:

[0028] κ m = span{r0, AP -1 r0,…,(AP -1 ) n-1 r0}

[0029] Where κ m represents an m-dimensional subspace, r0 represents the initial residual, and AP -1 represents the preconditioning matrix.

[0030] Furthermore, the iterative formula is as follows:

[0031] x m = x0 + P -1 V m y m

[0032] where x m represents an arbitrary vector, x0 represents the initial solution vector, V m represents the subspace basis, and y m represents the linear combination coefficient.

[0033] Furthermore, the process of optimizing the spatial tetrahedral mesh based on the quadratic nesting method includes:

[0034] Performing long-period edge expansion and Dirichlet boundary edge expansion on the spatial tetrahedron to construct a quadratic nested mesh.

[0035] Furthermore, the third-order backward Euler difference scheme is as follows:

[0036]

[0037] where M is the mass matrix, S is the stiffness matrix, J s is the field source, i, j represent matrix or vector elements, t represents time, and E j represents the electric field value on the edge, and n represents the time step after time discretization.

[0038] Furthermore, the process of reducing the computational error of iterative solution by using the mesh pre-refinement method includes:

[0039] Before the iterative solution, perform binary refinement on the field source and measurement points of the quadratic nested mesh to reduce redundant elements.

[0040] Furthermore, the process of implementing parallel computing using the domain decomposition method includes: adopting the domain decomposition method to divide the quadratic nested mesh into multiple regions, performing independent calculations on each region using MPI in parallel, after each time step calculation is completed, communicating between regions through MPI, restoring the local solution of each region to the global solution, and then performing the next iterative operation.

[0041] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0042] The present invention realizes the iterative solution of the TEM backward Euler finite element method based on the FGMRES algorithm, which not only ensures the computational efficiency but also greatly reduces the required memory. Based on unstructured tetrahedral meshes, it can better fit the terrain. Meanwhile, by optimizing the meshes and time steps, the computational efficiency is accelerated and the computational accuracy is improved, enabling the fast calculation of three-dimensional transient electromagnetic fields on a personal computer. The present invention verifies the correctness of our algorithm through cross-validation with analytical solutions, frequency-domain transformation solutions, and previous results. In addition, it is also compared with the direct solution method to analyze the advantages and disadvantages of the two methods, and the applicability of the algorithm is tested. The experimental results show that the proposed method has higher stability and lower requirements for mesh modeling and time step settings, providing a very convenient forward modeling means for subsequent inversion algorithms. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] By reading the following detailed description of the preferred embodiments, various other advantages and benefits will become clear to those of ordinary skill in the art. The drawings are only for the purpose of showing the preferred embodiments and are not considered to be a limitation of the present invention. In the drawings:

[0044] Figure 1 is a schematic flow chart of the method in the embodiment of the present invention;

[0045] Figure 2 is a schematic diagram of the spatial tetrahedral Nédélec H element and its electric field distribution in the embodiment of the present invention;

[0046] Figure 3 is a schematic diagram of mesh division in the embodiment of the present invention, where (a) is a traditional mesh and (b) is a secondary nested mesh;

[0047] Figure 4 is a schematic diagram of a layered medium model in the embodiment of the present invention;

[0048] Figure 5 is the calculation result of the transient electromagnetic method response in the embodiment of the present invention, where (a1)-(a4) are the induced electromotive force values (b1)-(b4) are the relative errors at 200m, (c1)-(c4) are the relative errors at 400m, and (d1)-(d4) are the relative errors at 1000m;

[0049] Figure 6 is a schematic diagram of a massive anomaly model in the embodiment of the present invention;

[0050] Figure 7 is a schematic diagram of the average relative difference of the transient electromagnetic responses of all measurement points in the embodiment of the present invention;

[0051] Figure 8The complex anomaly body model in the embodiment of the present invention;

[0052] Figure 9 The schematic diagram of the field source and measuring point distribution in the embodiment of the present invention;

[0053] Figure 10 The transient electromagnetic response calculation results in the embodiment of the present invention, where (a) is ground measurement, (b) is airborne measurement, (c) is the relative difference between iteration and direct solution, and (d) is the relative difference between considering the anomaly body and not considering the anomaly body. Specific implementation manners

[0054] Hereinafter, the exemplary embodiments of the present disclosure will be described in more detail with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that this disclosure will be more thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art. It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other. Hereinafter, the present invention will be described in detail with reference to the drawings and in conjunction with the embodiments.

[0055] This embodiment proposes a three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method, as Figure 1 shown, including:

[0056] Based on the vector finite element method, tetrahedral space discretization is performed, and the electric field is assigned to the six edges of the spatial tetrahedron;

[0057] The backward Euler method is used to discretize the time step, a backward Euler finite element equation set is constructed, and a linear system composed of a large sparse matrix is assembled based on this equation set;

[0058] Based on the finite element equation set, a linear system composed of a large sparse matrix is assembled;

[0059] Based on the FGMRES algorithm, the linear system is iteratively solved to obtain the electric field values corresponding to the edges of the spatial tetrahedron elements;

[0060] Based on the quadratic nesting method, the spatial tetrahedral mesh is optimized, and based on the third-order backward Euler difference format, the time step is optimized;

[0061] The mesh pre-refinement method is adopted to reduce the calculation error of the iterative solution, and the domain decomposition method is used to realize parallel calculation.

[0062] The detailed implementation process of this embodiment is as follows:

[0063] The control equation of the backward Euler finite element method:

[0064] The propagation of electromagnetic waves follows Maxwell's equations, which can be expressed in the time domain as follows:

[0065]

[0066] Where, E represents the electromagnetic intensity, H represents the magnetic field intensity, and J s Represents the field source, t represents time, ε is the dielectric constant, σ is the conductivity, and μ is the magnetic permeability. Combining the above equations and eliminating the magnetic field H, we can obtain the following damped wave equation:

[0067]

[0068] For transient electromagnetic method, the detection frequency is low enough to meet the quasi-static condition, and the displacement current term can be ignored. At this time, the above equation degenerates into the diffusion equation:

[0069]

[0070] Vector finite element method:

[0071] To ensure that the control equation xx satisfies the tangential continuity condition of the electric field, this embodiment uses the first-order Nédélec H (curl) unit of the Sobolev space tetrahedron for spatial discretization. The electric field is assigned to the six edges of the tetrahedron, such as Figure 2 As shown in the figure, the properties of the unit are given by the conductivity and magnetic permeability parameters. Using the first-order vector shape function and the tangential component value on the edge, the electric field value at any position in the unit can be expressed as:

[0072]

[0073] Where E represents the electric field value at any position, E j Represents the electric field value on the edge, N j represents the first-order interpolation basis function. Using Faraday's law of electromagnetic induction, the magnetic induction intensity can be expressed as:

[0074]

[0075] Where B represents the magnetic induction intensity and t represents time.

[0076] For any element in space, the following finite element equations are constructed using the Galerkin weighted residual method and the second-order backward Euler scheme:

[0077]

[0078] Where M is the mass matrix, S is the stiffness matrix, and J s is the field source, i, j represents the matrix or vector element (local edge number), the specific form is:

[0079]

[0080] Using the above formula, the edge numbers [i], [i,j] in all elements are assembled into the linear system matrix A according to the corresponding global edge numbers [Ndof i , [Ndof i ,Ndof j to obtain the linear system in matrix form:

[0081] Ax = b

[0082] where A = 3M + 2ΔtS, x = e n , b = M(4e n-1 - e n-2 ). For the assignment of the source term, the line source of the actual model can be discretized into several parts, and then each small segment of the wire is approximated as an electric dipole, and the field source is directly loaded on the edge (see Figure 2 ).

[0083] For the actual current waveform, the difference derivative is used to obtain the derivative of the current with respect to time, and finally it is directly substituted into the equation. For the simulation of the step response, in this embodiment, a long-period trapezoidal wave is used, and the switching function is used to adjust the rising edge and the falling edge to make the calculation relatively more stable.

[0084] Solution of the system of equations:

[0085] For the solution of large sparse systems of equations, direct solution and iterative solution methods can be used. For direct solution, in the control equation Ax = b, the left-hand side term only changes with the change of time. When the time step is constant, iterative solution can be directly carried out without repeating the decomposition of the left-hand side term. By using a piecewise constant increasing time step (Um), the calculation efficiency can be greatly improved. In the current finite element methods for calculating TEM, direct solution methods are all used. Although the calculation efficiency is relatively high, since the left-hand side term is a large sparse matrix, decomposing it will occupy a large amount of memory and require high computing equipment. Currently, many direct solvers have been developed according to different LU decomposition methods. Commonly used open-source solvers include MUMPS, PARDISO, SPOOLES, SuperLU, etc. In the subsequent comparative experiments, this embodiment uses the MUMPS open-source direct solver.

[0086] For the iterative solution method, although independent operations are performed for each time step, the time step can be flexibly changed, and the memory occupied is relatively small. The iterative solution has high requirements for the condition number of the equation. If the left-hand side term is not preprocessed, the iteration is difficult to converge. According to different iterative formats, different iterative solution methods are derived. Commonly used methods include stationary iterative methods such as Jacobi, Gauss-Seidel, and SOR, and subspace iterative methods such as CG, MINRES, GMRES, and BiCGStab. In the subsequent iterative calculation, the FGMRES iterative solver in the open-source computing library Hypre is used in this embodiment.

[0087] Iterative solution method based on right preconditioned generalized minimal residual method:

[0088] For the bi-curl equation system, the equation system is ill-conditioned, and the condition number of the matrix cond(A) >> 1. It is difficult to converge directly using the iterative solution method and the calculation accuracy is relatively low. After introducing preprocessing, this problem can be improved. To effectively reduce the condition number, a right preconditioner is added to the linear system equation to transform it into an equivalent equation system:

[0089] AP -1 Px = b

[0090] By setting a reasonable preconditioning matrix P, the condition number of the matrix (AP - 1) is much smaller than that of A, and the iterative solution can converge faster. Here, the auxiliary space Maxwell preconditioner (AMS) is used in this embodiment. This preconditioner has been proven to enable the equation to converge faster in the solution of the bi-curl equation. Based on this preconditioner, the calculation process of the FGMRES algorithm is as follows:

[0091] For the equation AP -1 Px = b, define an m-dimensional subspace

[0092] κ m = span{r0, AP -1 r0,..., (AP -1 ) n-1 r0}

[0093] where r0 represents the initial residual.

[0094] Let V m be a set of bases for the above space κ m , then any vector x can be expressed as x = V m y, where y is the linear representation coefficient. Using the iterative formula x m = x0 + P -1 Δx, Δx ∈ κ m , the problem is transformed into constructing the basis V m , and the linear representation coefficient y of Δx in this basism The calculation of which has an iterative formula at this time:

[0095] x m = x0 + P -1 V m y m

[0096] where x0 represents the initial solution vector.

[0097] Since it is very difficult to calculate P -1 , the inversion process in P -1 V m y m can be converted into solving the equation P -1 z i = v i (i = 1, …, m). At this time, using Z m = P -1 V m , the above formula can be converted into:

[0098] x m = x0 + Z m y m

[0099] At this time, the basis Z m = [z1, z2, …, z m is the extension of the subspace κ m . We use the Arnoldi process to construct an orthogonal basis. After obtaining the subspace bases V m , Z m , the parameters of the equation can be expressed as:

[0100] r = b - Ax m = r0 - AV m y = V m+1 (βe1 - H m+1 y)

[0101] where β = ‖r0‖2, The calculation of which is the minimization problem of the parameter r, that is

[0102] y m = arg min y ‖βe1 - H m+1,m y‖2

[0103] Using QR decomposition to solve the above formula, the linear representation coefficient y m is solved. Then, using the iterative formula x m = x0 + P -1 V m y m , the approximate solution x m can be solved. If the relative residual rm less than the specified initial value (e.g., 10 -6 ), then x = x m , otherwise let x0 = x m , and repeat the iteration.

[0104] Mesh generation and time step selection:

[0105] For the finite element method, the requirements for the spatial mesh are relatively high. However, blindly refining the mesh will increase the computational cost. For the calculation of the long-period transient electromagnetic response, the diffusion depth of the field can be expressed by the skin depth:

[0106]

[0107] In the formula, t obs represents the maximum observation time, μ0 represents the magnetic permeability of vacuum, and σ represents the conductivity of the half-space;

[0108] For a homogeneous half-space with a conductivity of 100 ohm-meters, the farthest diffusion range of the late-stage response up to 10 -1 seconds is about 4 kilometers, and it can reach 13 kilometers at 1 second. For layered media, if the calculation does not reach this range, the error of the late-stage calculation results will be relatively large. Usually, when using the Dirichlet boundary condition, the layer is extended outward in the extended boundary area, but this will result in relatively poor mesh quality in the extended boundary area of the layer (a relatively large mesh scale is required in the extended boundary area to reduce the spatial degrees of freedom), thereby increasing the calculation time. We reduce the number of meshes and improve the calculation accuracy by nesting the calculation area again. As Figure 3 shown, in subsequent calculations, except for the mesh test of the layered model, we use the ones that have been nested twice (the first time is for long-period extended boundary, and the second time is for Dirichlet boundary extended boundary).

[0109] For the selection of the time step, the algorithm we use is more flexible in the selection of the time step. In order to achieve high-precision calculation under the o(Δt 2 ) error, the step change rate still needs to be restricted to a certain extent (mainly depending on the time interval for which precise calculation is desired). We give the time step settings used in subsequent simulations. In fact, the time step can be further optimized to reduce the overall number of calculation time steps. When the time step changes greatly, we recommend using the third-order backward Euler format with o(Δt 3 ) accuracy (only storing one more matrix of the edge size, which has a relatively small impact on memory). The specific formula is as follows:

[0110]

[0111] where n represents the time step after time discretization.

[0112] Domain decomposition and mesh refinement:

[0113] For the iterative solution method, it is impossible to use a time-step format similar to the direct solution to accelerate the calculation, and the time consumption of direct calculation is relatively high. Considering the low memory occupancy characteristic of the iterative solution, in this embodiment, the domain decomposition method is used to divide the global mesh into multiple regions, and MPI parallelism is utilized to calculate separately on each region. After each time step calculation is completed, communication between regions is carried out through MPI to restore the local solution to the global solution, and then the next iterative operation is performed.

[0114] To ensure the calculation accuracy, it is usually necessary to accurately divide and refine the mesh, which requires high requirements for modeling. When facing too many measurement points, it will cause unit redundancy. In traditional frequency-domain electromagnetic algorithms, adaptive refinement is usually used to improve the calculation accuracy while solving quickly. However, for time-domain algorithms, it is difficult to define a suitable adaptive refinement formula. For this reason, this embodiment has conducted a large number of tests on the mesh and found that the units mainly affecting the accuracy are concentrated near the field source and measurement points, and only appropriate refinement is required for the underground space. Based on this, when actually performing calculations in this embodiment, on the basis of the original mesh, further binary refinement is performed on the field source and measurement points. Generally speaking, performing 3 to 5 times of refinement on the mesh directly before calculation can obtain good results, avoiding the multiple repeated calculations of optimizing the mesh by calculating the posterior error.

[0115] The present invention directly calculates the time-domain TEM response through the above steps, and this algorithm has the following advantages:

[0116] Based on the AMS preconditioner, the matrix condition number can be effectively reduced, and the stability of the Krylov method can be increased; through the iterative solution algorithm, the calculation memory can be effectively reduced, the requirements for the calculation device are relatively low, and at the same time, forward modeling can be performed on large-scale geological models. Based on the domain decomposition parallel algorithm, on different calculation devices, the calculation efficiency can be improved or the memory requirement can be reduced by modifying the number of processes.

[0117] Verification based on the layered medium model:

[0118] In the first verification example, this embodiment calculates the numerical solution and the one-dimensional semi-analytical solution of the layered model (as shown in Fig. 1) to verify the accuracy of the algorithm. Since the one-dimensional semi-analytical solution is directly calculated from the frequency-domain analytical solution and converted to the time domain, it has very high accuracy. Therefore, we define the forward relative error as:

[0119]

[0120] In the formula, e represents the relative error, represents the one-dimensional semi-analytical solution, Represents the numerical solution.

[0121] Consider a four-layered model as Figure 4 shown. The corresponding resistivities from top to bottom are (500, 100, 50, 10) ohm-meters respectively, the layer thicknesses are (200, 300, 500) meters respectively, a wire source with a length of 500 meters is located at (0, -250 to 250, 0) meters, the offsets of three receiving points arranged along the x-axis are (200, 400, 1200) meters respectively, the total degrees of freedom of the model are 382168, and the time-domain response range is 1e -6 ~1e -1 s.

[0122] We calculated the TEM responses at the time steps of 310 and 120 with the backward Euler orders of 2 and 3 respectively. Meanwhile, we also calculated the results based on the MUMPS direct solver, which are at the time steps of 120 with 6 factorizations and at the time steps of 310 with 10 factorizations respectively. For the time sampling of the direct solution, we first determine the time nodes at logarithmically equidistant time steps, and then insert the iterative time steps at equal intervals between adjacent time nodes. Finally, we also calculated the TEM response under the conventional extended-edge grid, and the degrees of freedom of this grid are 732567.

[0123] Figure 5 Shows the calculation results and relative errors corresponding to seven different calculation schemes. The black solid line represents the one-dimensional semi-analytical solution, and the blue and red solid lines represent the numerical solutions. Generally speaking, all six schemes can accurately calculate the transient electromagnetic response results, and the relative errors are all below five percent, meeting the accuracy requirements.

[0124] Table 1

[0125]

[0126] Table 1 shows the computational resources consumed corresponding to seven schemes. Comparing the time discretization schemes of 120 and 310 steps, it can be seen that for iterative solution, only dozens of time steps are needed for accurate calculation. Although further densifying the time steps can slightly improve the calculation accuracy, it requires more computational time. Comparing the calculation results of the second- and third-order backward Euler schemes, since the second-order scheme can already perform accurate calculations, the improvement of the third-order scheme for calculation error is relatively limited. However, the calculation times and memory occupancies of the second and third orders are almost the same. To ensure high-precision simulation of complex models, the third-order time accuracy scheme is still recommended. Comparing the results of the regular grid and the quadratic extended-edge grid, due to more densification of the main computational region, the quadratic extended-edge grid has higher calculation accuracy. At the same time, due to lower degrees of freedom, its calculation speed is also accelerated by about 1.6 times. Comparing the direct and iterative solvers, for the parallel mumps solver, by default, all cores are used to accelerate the calculation speed. However, for a personal computer, due to the limitation of the CPU, the domain decomposition algorithm cannot fully reflect the advantage of multi-core parallel computing. When the number of cores reaches 6 cores, the calculation efficiency no longer improves. Nevertheless, in this embodiment, the computational resource consumptions in the same environment are still compared to prove the advantage of the algorithm. When switching the step size, the result error of the direct solution method is relatively large. When the time step of the direct solution method is densified to a certain extent, its result is very similar to that of the iterative solution method. In this example, the direct solution method only needs 10 decompositions to obtain relatively accurate results. The iterative method has lower requirements for time steps and can perform calculations with fewer time steps. Under the same accuracy conditions, the calculation efficiency is increased by about 1 time compared to the direct solution method, and the memory usage is reduced by about 1.7 times.

[0127] Verification is carried out based on the block anomaly model:

[0128] In the second verification example, to further verify the applicability of the iterative solution method and its corresponding grid construction and time step scheme, a block anomaly model is used. The model is composed of three block anomalies under a homogeneous half-space, as Figure 5 shown.

[0129] Table 2 gives the vertex information corresponding to the lower left and upper right corners of the three block anomalies. A wire source with a length of 500 meters is located at (0, -250 to 250, 0) meters. Three survey lines are set on the surface (z = 0 meters) on its right side, and the corresponding x coordinates are (-1000, 0, 1000) meters. Six measurement points are set on each survey line, and the corresponding y coordinates are (500, 1000, 1500, 2000, 2500, 3000) meters.

[0130] Table 2

[0131]

[0132] There is not much optimization for the establishment of the grid. Only the measuring points and the positions of the field sources are further encrypted and optimized in the calculation, and the coefficient construction of the entire grid is relatively simple. By controlling the growth ratio of the grid and the maximum element size, the accuracy of the calculation is ensured. Table 3 gives some construction parameters of the grid. In fact, all calculation results use a similar growth ratio, and most of them modify the maximum element volume in different spatial domains, with fewer changes in other parameters.

[0133] Table 3

[0134]

[0135] In this embodiment, the TEM responses for 120 steps under the 2nd and 3rd order backward Euler formats are calculated, and the results for 310 steps under the 2nd order format are calculated for comparison. To verify the correctness of the calculation results, they are compared with the calculation results of the same model in the open-source program custem. To calculate the relative difference between two non-exact numerical solutions, the relative difference adopted is defined as:

[0136]

[0137] where e represents the relative error, represents the result of the open-source software, represents the numerical solution.

[0138] Figure 6 They are the calculation results and relative errors of five calculation methods at the measuring point (1000, -1000, 0) meters. It can be seen that the calculation accuracy is lower than 5%, meeting the calculation error requirements. By comparing the results of the second-order and third-order formats for 120 steps, it can be seen that the third-order format has higher calculation accuracy for complex models, especially in the late-time channels. By comparing the second-order calculation results with different time steps, it can be seen that the results after encrypting the time have higher accuracy, and encrypting the time step can improve the calculation accuracy to a certain extent. By comparing the calculation results with different step lengths and orders, there is no significant difference in the calculation errors between the two modes. This embodiment also gives the results of direct solution. The relative error for 120 steps is greater than 5%, and 310 steps are required to meet the error requirements. Using the third-order time format can reduce the requirements for the time step and has stronger applicability.

[0139] Figure 7The average errors of all measurement points under the second-order (a) and third-order (b) calculation methods are shown. It can be seen that compared with the second-order format, the calculation error of the third-order format at the measurement points close to the field source is significantly reduced. For long wire sources, strong equation singularities will occur at positions close to the field source and where resistivity changes are complex. The spatial error can be reduced by grid refinement and using high-order elements, and the time error can be reduced by densifying the time steps and using high-order time discretization formats.

[0140] Table 4

[0141]

[0142] Table 4 (computational efficiency and memory occupancy of the massive anomaly model) shows the computational resource consumption of five calculation modes. We also list the consumption when using open-source software to calculate the same model. Compared with the second-order format, the time and memory usage for the accuracy of the third-order format are basically the same. Compared with the direct solver, the iterative solver reduces the memory consumption required for calculation, and due to its weaker dependence on the time step, longer time steps can be used to improve the computational efficiency. And because the scale of this model is relatively larger, the relative improvement in computational efficiency is greater.

[0143] Complex large-scale model:

[0144] From the first and second examples, it can be seen that the optimization of the spatial grid modeling and time step selection in the method described in this embodiment can effectively improve the computational efficiency. For transient electromagnetic methods, the actual underground medium model is very complex, and accurate modeling of complex geological structures is required. However, the calculation of TEM responses for large and complex models has a very high memory occupancy and relatively low computational efficiency. To further verify the scalability of the algorithm, in this embodiment, a three-dimensional complex geological body model is established according to the actual environment, such as Figure 8 shown. The main body of the model is four inclined layers considering the actual geological structure. There is an upward low-resistance fracture zone and a low-resistance abnormal ore body within the layers. The corresponding resistivity values are marked in the figure. The model is constructed based on the actual exploration data of the Dongguashan area in Tongling. Here, the resistivity is only used for forward modeling tests and is not the real resistivity.

[0145] A wire source with a length of 1 km is arranged along the x-axis direction. The source is divided into 50 parts and discretely assigned according to the actual terrain data. A total of 41 survey lines are arranged along the y-axis direction. The x-coordinate ranges from 500 m to 2500 m, with a line spacing of 50 m. Each survey line has 21 measurement points. The y-coordinate ranges from 1000 m to 3500 m, with a line-point spacing of 125 m, totaling 861 measurement points, as Figure 9As shown. The calculation results collected by both ground and semi-aerial methods were tested separately to highlight the influence of terrain on the transient electromagnetic method. Among them, the ground measurement points were distributed along the actual terrain, and the height of the aerial measurement points from the horizontal plane was 300 meters.

[0146] Due to the high requirement of the complex model for the time step, the time step was encrypted during the forward modeling in this experiment. Here, the computing CPUs of 16 cores and 32 cores are Intel(R) Xeon(R) Gold 6248R, and the others are i5-13500HX. The iterative solution time step used 310 steps mentioned above. For the direct solution, the relative difference of 310 steps was relatively large, so 20 decompositions were adopted and the time step was encrypted to 1312 steps. At the same time, the symmetric mean absolute percentage error of the two methods was calculated to verify the stability of the algorithm.

[0147] Figure 10 The transient electromagnetic responses of the survey lines at 0.164 s for the ground (a) and the air (b) are shown. By comparing the measurement results of different observation methods, it can be clearly seen that the terrain has a greater impact on the data. The measurement results of the semi-aerial measurement are smoother, which is beneficial to data processing, analysis and inversion. We calculated the relative errors of the two schemes, as Figure 10 (c) shows, less than 5%, which also verifies the correctness of the algorithm in this embodiment. In addition, the relative influence of the abnormal body on the field value under the semi-aerial device was calculated, and the results are shown in Figure 10 (d). It can be seen that for the low-resistivity ore body of 200 - 400 meters, its influence on the field value reaches 10%, which can be observed in actual exploration.

[0148] Table 5

[0149]

[0150] Table 5 shows the time consumption and memory occupancy of the iterative solution and direct solution schemes under different numbers of parallel processes. For the two solution schemes, under the condition of sufficient memory and CPU cores, the iterative solution scheme of this embodiment can greatly improve the calculation efficiency by enabling multi-core parallelism. At the same time, under the same computing device, the iterative solution has a faster calculation speed and lower memory occupancy, and the program has stronger scalability.

[0151] The present invention proposes a subspace iteration solution method using the Generalized Minimum Residual method (Fgmres) and utilizes domain decomposition for parallel acceleration. Specifically, for the calculation of long-period transient electromagnetic responses, the grid is optimized, a set of modeling parameters with strong adaptability is selected, and the grid degrees of freedom are reduced without affecting the accuracy. In addition, the present invention replaces the commonly used second-order format with a third-order backward Euler difference format, conducts forward modeling tests with different time steps, gives a relatively reasonable time step, and improves the calculation efficiency. By comparing the layered model with the semi-analytical solution and the three-dimensional model with the results of predecessors, the effectiveness of the algorithm is verified. The comparison results with other open-source codes and the numerical experiments on large-scale models show that this method can achieve high accuracy and efficiency for both simple and complex models, and at the same time occupies relatively less memory.

[0152] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the above embodiments, those of ordinary skill in the art should understand that: modifications or equivalent replacements can still be made to the specific embodiments of the present invention, and any modifications or equivalent replacements that do not depart from the spirit and scope of the present invention should be covered within the protection scope of the claims of the present invention.

Claims

1. A three-dimensional finite element transient electromagnetic forward modeling calculation method based on an iterative solution method, characterized in that, Including: Based on the vector finite element method, tetrahedral space discretization is carried out, and the electric field is assigned to the six edges of the space tetrahedron; The backward Euler method is used to discretize the time step, a backward Euler finite element equation set is constructed, and a linear system composed of a large sparse matrix is assembled based on the backward Euler finite element equation set; Based on the FGMRES algorithm, iterative solution of the linear system is carried out to obtain the electric field values corresponding to the edges of the space tetrahedron elements; Based on the quadratic nesting method, the space tetrahedron mesh is optimized, and the time step is optimized based on the third-order backward Euler difference format; The mesh pre-refinement method is used to reduce the computational error of iterative solution, and the domain decomposition method is used to achieve parallel computing.

2. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that, The finite element equation set is as follows: where M is the mass matrix, S is the stiffness matrix, J s is the field source, i, j represent the matrix or vector elements, t represents time, E j represents the electric field value on the edge, and n represents the time step after time discretization.

3. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that, The linear system matrix is as follows: Ax = b where A represents the linear system matrix, A = 3M + 2ΔtS, x = e n , b = M(4e n-1 - e n-2 ), M is the mass matrix, S is the stiffness matrix, and e n represents the electric field after time discretization.

4. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that The process of iterative solution of the linear system matrix based on the FGMRES algorithm includes: Based on the auxiliary space Maxwell preconditioner, a preconditioned matrix is constructed for the linear system matrix to reduce the matrix condition number, and the linear system matrix is converted into an equivalent equation set; An m-dimensional subspace is defined for the equivalent equation set, and a vector at any position within the space tetrahedron element is represented based on a set of bases in the m-dimensional subspace; An iterative formula is constructed based on the linear representation coefficients, the m-dimensional subspace is expanded, a subspace basis is constructed, and the linear representation coefficients are solved based on the subspace basis; Based on the linear representation coefficients, the iterative formula is calculated to obtain the vector representation at any position within the space tetrahedron element.

5. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 4, characterized in that The m-dimensional subspace is as follows: κ m = span{r0, AP -1 r0, …, (AP -1 ) n-1 r0} where κ m represents an m-dimensional subspace, r0 represents the initial residual, and AP -1 represents the preconditioning matrix.

6. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 4, characterized in that The iterative formula is as follows: x m = x0 + P -1 V m y m where x m represents an arbitrary vector, x0 represents the initial solution vector, and V m represents the subspace basis, and y m represents the linear representation coefficient.

7. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that The process of optimizing the space tetrahedron mesh based on the quadratic nesting method includes: Long-period edge expansion and Dirichlet boundary expansion are carried out on the space tetrahedron to construct a quadratic nested mesh.

8. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that The third-order backward Euler difference format is as follows: Among them, M is the mass matrix, S is the stiffness matrix, J s is the field source, i, j represent matrix or vector elements, t represents time, E j represents the electric field value on the edge, and n represents the time step after time discretization.

9. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that The process of using the mesh pre-refinement method to reduce the computational error of iterative solution includes: Before the iterative solution, the field sources and measurement points of the quadratic nested mesh are bisected and encrypted to reduce redundant elements.

10. The three-dimensional finite element transient electromagnetic forward calculation method based on the iterative solution method according to claim 1, characterized in that, The process of using the domain decomposition method to achieve parallel computing includes: adopting the domain decomposition method to divide the quadratic nested mesh into multiple regions, performing independent calculations on each region using MPI parallelism, after each time step calculation is completed, communicating between regions through MPI, restoring the local solutions of each region to the global solution, and then performing the next iterative operation.

Images