Transient electromagnetic rapid forward modeling method based on integral restart algorithm

Abstract

The invention discloses a transient electromagnetic fast forward modeling method based on an integral restart algorithm, which adopts a vector finite element method discrete control equation of a tetrahedral mesh, and writes a secondary field at any moment after a current is turned off into a form of multiplying a matrix exponential function f (A) by a vector b. Firstly, calculation of a matrix exponential function f (A) is converted into an integral form by using a Gauss-Legendre numerical integration technology, then approximation is performed on the integral by using an Arnoldi algorithm, and finally rapid forward modeling is realized. According to the method, the convergence speed is high, the memory occupation is low, and the instability problem caused by error accumulation in the Arnoldi algorithm restarting process can be avoided.

Description

Technical Field

[0001] The present invention belongs to the field of transient electromagnetic detection, and more specifically, it is a fast forward modeling method for transient electromagnetic based on an integral restart algorithm. Background Art

[0002] The transient electromagnetic method uses grounded electrodes or ungrounded loops to transmit bipolar pulsed currents into the ground. The induced eddy currents in the underground medium generate a secondary field that varies with time. During the intermittent period of the primary field, a receiving coil is used to measure the magnetic field signal. By extracting and analyzing the secondary field signal, the purpose of detecting underground geological bodies can be achieved. Currently, as one of the main non-seismic methods, it is widely used in near-surface resource, environmental, and engineering explorations. Forward modeling is crucial for the analysis of the response characteristics of typical abnormal bodies and the inversion and interpretation of transient electromagnetic data. However, since the order of the coefficient matrix A generated by three-dimensional transient electromagnetic forward modeling often reaches hundreds of thousands or even millions, it is very difficult to use existing numerical methods to have a satisfactory calculation efficiency while occupying a small amount of memory. That is, in the current processing method, if better calculation efficiency is required, a large amount of memory needs to be occupied; if only a small amount of memory is occupied, the calculation efficiency is very low.

[0003] Therefore, how to provide a new method that can maintain a small amount of memory occupied and still have a high calculation efficiency during the transient electromagnetic forward modeling process is the research direction required by the present invention. Summary of the Invention

[0004] In view of the problems existing in the above-mentioned prior art, the present invention provides a fast forward modeling method for transient electromagnetic based on an integral restart algorithm, which can maintain a small amount of memory occupied and still have a high calculation efficiency during the transient electromagnetic forward modeling process.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is: a fast forward modeling method for transient electromagnetic based on an integral restart algorithm, and the specific steps are as follows:

[0006] Step 1: Model parameterization: Establish a geophysical geoelectric model according to known geological data;

[0007] Step 2: Spatial discretization processing: Use the vector finite element method to discretize the governing equation to obtain the matrix exponential function f(A)b required after spatial discretization processing;

[0008] Step 3: Determine the integral expression of the error function: First, use an interpolation polynomial to process the function f in Step 2, and then determine the Arnoldi approximation f m , so as to obtain the error between the Arnoldi approximation f m and f(A)b as the error matrix function e m (A) multiplied by a vector;

[0009] Step 4. Arnoldi restart calculation of the approximation value: Use the Arnoldi restart algorithm to process the error matrix function e obtained in Step 3 m (A) First, perform the first processing to obtain the Arnoldi approximation value after the first restart loop Then, use the Arnoldi approximation value Continue to perform subsequent multiple restart loop processes to obtain the Arnoldi approximation values for each subsequent restart process

[0010] Step 5. Data forward modeling: According to the Arnoldi approximation values for each restart loop calculated in Step 4 Calculate the corresponding matrix exponential function f(A)b, and finally use the matrix exponential function f(A)b and perform forward modeling simulation on the geoelectric model in Step 1 through GPU parallel processing.

[0011] Furthermore, the specific content of Step 2 is as follows:

[0012] The time-domain Maxwell equation ignoring the displacement current is as follows

[0013]

[0014] In the above formula, E(r,t) is the electric field strength (V / m), B(r,t) is the magnetic induction intensity (T), t is the time (s), μ = μ 0 is the vacuum permeability (H / m), σ(r) is the conductivity of the medium in the calculation area (S / m), s(r,t) is the external source term (A / m 2 );

[0015] Since the downward step wave is widely used in transient electromagnetic detection, the present invention uses the downward step wave as the emission waveform; the external current source generated by the transient electromagnetic instrument is turned off at t = 0;

[0016] s(r,t) = q(r)H(-t) (2)

[0017] Here, H is the unit Heaviside function:

[0018]

[0019] The vector field q(r) is the discrete current density;

[0020] For the downward step waveform, at time t = 0, only a steady magnetic field exists in the whole space, that is, the initial condition;

[0021] B(r,t = 0) = B 0 (4)

[0022] Initial value magnetic induction intensity B 0 Obtained through the analytical expression;

[0023] Assume that the boundary is located in a sufficiently far region, apply the first-kind boundary condition to the boundary in the spatial domain, and assume that the electromagnetic field can decay to zero, that is

[0024] B×n = 0; n×E = 0 (5)

[0025] Discretize the governing equations using the vector finite element method on the tetrahedral mesh and eliminate the electric field, then Equation (1) can be rewritten as a system of ordinary differential equations with an initial value problem;

[0026]

[0027] In Equation (6), the coefficient matrix is a large sparse asymmetric matrix, C is the curl operator matrix, M σ and M μ are matrices containing conductivity and mesh information respectively; the analytical solution of the time-domain electromagnetic equation (6) is expressed using the matrix exponential function, specifically as follows:

[0028] B(t) = exp(-tA)B 0 (7)

[0029] And let f(A) = exp(-tA), b = B 0 , to obtain the matrix exponential function f(A)b.

[0030] Furthermore, the specific content of the third step is as follows:

[0031] First, represent the function f in the form of a Cauchy-type integral, and use the interpolation polynomial to estimate the value of the function f at the given nodes; let Ω be a region in the real number field and let be an analytical function, which has the following integral form, the path and the function

[0032]

[0033] The interpolation polynomial p of f with interpolation nodes is: m-1 is:

[0034]

[0035] Here w m (z) = (z - θ 1 )…(z - θ m );

[0036] Next, the error integral expression of the Arnoldi approximation \(f(A)b\) is given: Since the matrix \(\text{spec}(A)\in\Omega\), and the vector is known, then using \(F\) m to represent the \(m\)-th Arnoldi approximation of \(f(A)b\);

[0037] Since \(p\) m-1 is the interpolation polynomial of \(f\) with interpolation nodes \(\theta\) 1 , \(\theta\) 2 , …, \(\theta\) m , subtracting \(p\) m-1 from the expression (11), we get:

[0038]

[0039] Replacing \(z\) in formula (10) with the matrix \(A\), and then multiplying by the vector \(b\); since \(p\) m-1 (A)b = f m , we get:

[0040]

[0041] Here

[0042] From the above formula, the error between the Arnoldi approximation \(f\) m and \(f(A)b\) is the product of an error matrix function \(e\) m (A) and a vector.

[0043] Furthermore, the specific content of the fourth step is as follows:

[0044] For the integration node \(t\) i \(\in\Gamma\) and the weight factor \(\omega\) i , any numerical integration formula of \(e\) m (A) in formula (11) has the following form

[0045]

[0046] Let the integration nodes and weight factors in formula (12) be fixed in all restart cycles; in the first restart loop, i.e., \(w\) m (t i ) = 1, then the first Arnoldi approximation value can be obtained:

[0047]

[0048] In subsequent restart processes, i.e., \(k\geq2\), the approximation value of each restart loop can be calculated by using the error expression (12)

[0049]

[0050] Finally, the following expression is obtained

[0051]

[0052] Furthermore, the specific content of Step 5 is as follows:

[0053] The Arnoldi approximation value of each restart cycle calculated according to Step 4 Thereby, the corresponding matrix exponential function f(A)b is calculated. Since f(A)=exp(-tA) and b = B 0 , since the Arnoldi algorithm is essentially equivalent to an iterative algorithm with adaptive time step selection, in order to accelerate the calculation speed, formula (7) is rewritten as

[0054] B(t i ) = exp(-h i A)B i-1 (18)

[0055] Here, t i = t i-1 + h i , i = 1, 2, …, nt, h i is the time step, and nt is the number of discrete time channels;

[0056] Then, the forward modeling simulation of the geoelectric model in Step 1 is carried out by using the rewritten formula (18) and through the GPU parallel processing method; during the GPU parallel processing process, the specific processing steps are as follows: First, copy the matrix data from the main memory to the GPU memory, guide the calculation process through the CPU, then quickly execute intensive numerical operations on each core of the GPU, and finally return the calculation results to the main memory to complete the forward modeling process.

[0057] Compared with the prior art, the present invention discretizes the control equation by using the vector finite element method of tetrahedral meshes, and writes the secondary field at any time after the turn-off current in the form of the multiplication of a matrix exponential function f(A) and a vector b. However, since the order of the coefficient matrix A generated by the transient electromagnetic forward modeling often reaches hundreds of thousands or even millions, it is almost impossible to use traditional numerical methods to still have satisfactory computational efficiency while occupying less memory. The present invention utilizes the integral expression of the iterative error of the Arnoldi method and successfully develops an efficient restart algorithm. First, the Gauss-Legendre numerical integration technique is used to transform the calculation of the matrix exponential function f(A) into an integral form, and then the Arnoldi algorithm is used to approximate this integral, and finally the forward modeling calculation is realized. The present invention not only has a fast convergence rate and low memory occupancy, but also can avoid the instability problem caused by error accumulation during the restart process of the Arnoldi algorithm; through numerical simulation verification, during the transient electromagnetic forward modeling process of the present invention, it can have high computational efficiency while maintaining small memory occupancy. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 is the 3D complex geoelectric model set in the numerical simulation verification.

[0059] Figure 2 is the calculation result of the electromagnetic field in the numerical simulation verification.

[0060] In the figure: (a) magnetic induction intensity containing a low-resistivity anomaly; (b) magnetic induction intensity without a low-resistivity anomaly; (c) relative anomaly. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0061] The present invention will be further described below.

[0062] Embodiment: The specific steps are as follows:

[0063] Step 1, model parameterization: According to the known geological data, establish a geophysical geoelectric model;

[0064] Step 2, spatial discretization processing: Use the vector finite element method to discretize the control equation to obtain the matrix exponential function f(A)b required after spatial discretization processing, specifically:

[0065] The time-domain Maxwell equation ignoring the displacement current is as follows

[0066]

[0067] In the above formula, E(r,t) is the electric field strength (V / m), B(r,t) is the magnetic induction intensity (T), t is the time (s), μ = μ 0is the vacuum permeability (H / m), σ(r) is the conductivity of the medium in the calculation region (S / m), and s(r,t) is the external source term (A / m 2 );

[0068] Since the downward step wave is widely used in transient electromagnetic detection, the present invention uses the downward step wave as the emission waveform; the external current source generated by the transient electromagnetic instrument is turned off at t = 0;

[0069] s(r,t) = q(r)H(-t) (2)

[0070] Here, H is the unit Heaviside function:

[0071]

[0072] The vector field q(r) is the discrete current density;

[0073] For the downward step waveform, at time t = 0, only a steady magnetic field exists in the entire space, that is, the initial condition;

[0074] B(r,t = 0) = B 0 (4)

[0075] The initial value of the magnetic induction intensity B 0 is obtained through an analytical expression;

[0076] Assume that the boundary is located in a sufficiently far region. Apply the first-kind boundary condition to the boundary in the spatial domain and assume that the electromagnetic field can decay to zero, that is

[0077] B×n = 0; n×E = 0 (5)

[0078] Discretize the control equation using the vector finite element method on a tetrahedral mesh and eliminate the electric field, then Equation (1) can be rewritten as an ordinary differential equation system with an initial value problem;

[0079]

[0080] In Equation (6), the coefficient matrix is a large sparse asymmetric matrix, C is the curl operator matrix, M σ and M μ are matrices containing conductivity and grid information respectively; the analytical solution of the time-domain electromagnetic equation (6) is expressed using the matrix exponential function, specifically as follows:

[0081] B(t) = exp(-tA)B 0 (7)

[0082] Let f(A) = exp(-tA), b = B 0 , and obtain the matrix exponential function f(A)b.

[0083] Step 3. Determine the integral expression of the error function: First, use the interpolation polynomial to process the function f in Step 2, and then determine the Arnoldi approximation f m , so as to obtain the Arnoldi approximation f m The error between f(A)b is the error matrix function e m (A) and the product of a vector, specifically:

[0084] First, express the function f in the form of a Cauchy-type integral, and use the interpolation polynomial to estimate the value of the function f at the given nodes; let Ω be the real number field in a region, and let be an analytic function, which has the following integral form, the path and the function

[0085]

[0086] has interpolation nodes The interpolation polynomial p of f m-1 is:

[0087]

[0088] where w m (z) = (z - θ 1 )…(z - θ m );

[0089] Next, give the error integral expression of the Arnoldi approximation f(A)b: Since the matrix spec(A) ∈ Ω, and the vector is known, then use f m to represent the m-th Arnoldi approximation of f(A)b;

[0090] Since p m-1 is the interpolation polynomial of f with interpolation nodes θ 1 , θ 2 ,…, θ m , use the expression (11) minus p m-1 , and get:

[0091]

[0092] Replace z in formula (10) with the matrix A, and then multiply by the vector b; since p m-1 (A)b = f m , get:

[0093]

[0094] Here

[0095] From the above formula, the Arnoldi approximation f m has an error of an error matrix function e with respect to f(A)b m which is the product of (A) and a vector.

[0096] Step 4. Arnoldi restart calculation of the approximation: Use the Arnoldi restart algorithm to first process the error matrix function e m (A) obtained in Step 3 to obtain the Arnoldi approximation after the first restart loop Then use the Arnoldi approximation to continue with subsequent multiple restart loop processes, thereby obtaining the Arnoldi approximations for each subsequent restart process Specifically:

[0097] For the integration nodes t i ∈Γ and the weight factor ω i , any numerical integration formula for e m (A) in formula (11) has the following form

[0098]

[0099] Let the integration nodes and weight factors in formula (12) be fixed in all restart cycles; in the first restart loop, i.e., w m (t i ) = 1, then the first Arnoldi approximation can be obtained:

[0100]

[0101] In subsequent restart processes, i.e., k ≥ 2, the approximation for each restart loop can be calculated by using the error expression (12)

[0102]

[0103] Finally, the following expression is obtained

[0104]

[0105] In this embodiment, the integration nodes are initialized and If the norm of the difference between the obtained integration approximations is greater than the specified error tolerance tol, then the number of integration nodes is further increased to Multiply by a factor until the desired numerical accuracy is achieved; if the number of integration nodes does not increase in a given restart cycle, start the next restart cycle with a reduced number of integration nodes.

[0106] Step Five, Forward Modeling of Data: Based on the Arnoldi approximation for each restart cycle calculated in Step Four calculate the corresponding matrix exponential function f(A)b, and finally use the matrix exponential function f(A)b and perform forward modeling simulation on the geoelectric model in Step One through GPU parallel processing, specifically as follows:

[0107] Based on the Arnoldi approximation for each restart cycle calculated in Step Four calculate the corresponding matrix exponential function f(A)b. Since f(A)=exp(-tA) and b = B 0 , because the Arnoldi algorithm is essentially equivalent to an iterative algorithm for adaptively selecting time steps, to speed up the calculation, rewrite formula (7) as

[0108] B(t i ) = exp(-h i A)B i-1 (18)

[0109] Here, t i = t i-1 + h i , i = 1, 2,..., nt, h i is the time step, and nt is the number of discrete time channels;

[0110] Then use the rewritten formula (18) and perform forward modeling simulation on the geoelectric model in Step One through GPU parallel processing; during the GPU parallel processing, the specific processing steps are as follows: First, copy the matrix data from the main memory to the GPU memory, guide the calculation process through the CPU, then quickly execute intensive numerical operations on each core of the GPU, and finally return the calculation results to the main memory to complete the forward modeling process.

[0111] Numerical Simulation Verification:

[0112] Set up a large-scale 3D complex geoelectric model to verify the method of this embodiment. The model schematic diagram and the resistivity of the strata are as Figure 1 shown. Set the resistivity of the air layer to 10 6Ω·m, the resistivity of the homogeneous half-space formation is 500 Ω·m. The surface is Quaternary sediments with a resistivity of 300 Ω·m; the resistivity of the shallow geothermal reservoir is 100 Ω·m, and the resistivities of the deep geothermal reservoirs are 50 Ω·m and 10 Ω·m respectively, and the resistivity of the geothermal fluid is 8 Ω·m. The lengths of these geological anomalies in the y direction are all 1 km. The length of the grounding wire source used in this paper is 6 km, and the coordinates of the two endpoints are (-3 km, -100 m, 0 m) and (3 km, -100 m, 0 m) respectively. The receiving survey line is parallel to the emission source, the offset is 200 m, the point spacing is 20 m, and the coordinates of the end of the survey line are (-2.5 km, 100 m, 0 m) and (2.5 km, 100 m, 0 m) respectively. The number of unknowns to be solved in this complex geoelectric model is as high as 6.86 million. The average calculation time of the existing method is 398.7 s, and the average memory consumption is 8.3 GB, while the calculation time using the present invention is only 32.4 s, and the memory consumption is only 3.1 GB. This shows that the present invention has high calculation efficiency on the premise of consuming less memory.

[0113] Figure 2 They are the electromagnetic responses measured on the surface with and without anomalies and the relative error between the two. From Figure 2 c, it can be seen that the forward simulation method in this paper can effectively reflect the electromagnetic anomaly characteristics of geological bodies.

[0114] The above are only the preferred embodiments of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. A transient electromagnetic fast forward modeling method based on integral restart algorithm, characterized in that: The specific steps are: Step 1: Model parameterization: Establish a geophysical and geoelectrical model based on geological data; Step 2: spatial discretization: use the vector finite element method to discretize the control equation and obtain the matrix exponential function f(A)b required after spatial discretization; Step 3: Determine the integral expression of the error function: First, use the interpolation polynomial to process the function f in step 2, and then determine the Arnoldi approximation f m , thus obtaining the Arnoldi approximation f m The error with f(A)b is the error matrix function e m (A) Multiplication with a vector; Step 4: Arnoldi restart calculation approximation: Use the Arnoldi restart algorithm to calculate the error matrix function e obtained in step 3 m (A) Perform the first processing first to obtain the Arnoldi approximation after the first restart cycle Then, using the Arnoldi approximation Continue to perform subsequent multiple restart cycles to obtain the Arnoldi approximation for each subsequent restart process Step 5: Data forward modeling: Arnoldi approximation for each restart cycle calculated in step 4 The corresponding matrix exponential function f(A)b is thus calculated, and finally the matrix exponential function f(A)b is used to perform forward simulation on the geoelectric model of step one through GPU parallel processing.

2. The transient electromagnetic fast forward modeling method based on the integral restart algorithm according to claim 1 is characterized in that: The step 2 is specifically as follows: Ignoring the displacement current, the time domain Maxwell equation is as follows In the above formula, E(r, t) is the electric field intensity, B(r, t) is the magnetic induction intensity, t is the time, μ = μ0 is the vacuum magnetic permeability, σ(r) is the conductivity of the medium in the calculation area, and s(r, t) is the external source term; Since the down-step wave is widely used in transient electromagnetic detection, the present invention adopts the down-step wave as the transmitting waveform; the external current source generated by the transient electromagnetic instrument is turned off at t=0; s(r,t)=q(r)H(-t) (2) Here H is the unit Heaviside function: The vector field q(r) is the discrete current density; For the lower step waveform, at time t = 0, there is only a steady magnetic field in the entire space, which is the initial condition; B(r,t=0)=B0 (4) The initial value of magnetic induction intensity B0 is obtained by analytical expression; Assuming that the boundary is located in a sufficiently far region, the first-type boundary condition is applied to the boundary in the space domain, and it is assumed that the electromagnetic field can decay to zero, that is, B×n=0; n×E=0 (5) By using the vector finite element method on a tetrahedral mesh to discretize the governing equations and eliminate the electric field, equation (1) can be rewritten as a system of ordinary differential equations with an initial value problem; In equation (6), the coefficient matrix is a large sparse nonsymmetric matrix, C is the curl operator matrix, M σ and M μ are matrices containing conductivity and grid information respectively; the analytical solution of the time domain electromagnetic equation (6) is expressed by a matrix exponential function, as follows: B(t)=exp(-tA)B0 (7) And let f(A)=exp(-tA), b=B0, and obtain the matrix exponential function f(A)b.

3. The transient electromagnetic fast forward modeling method based on the integral restart algorithm according to claim 2 is characterized in that: The step three is specifically as follows: First, the function f is expressed as a Cauchy-type integral, and the interpolation polynomial is used to estimate the value of the function f at a given node; let Ω be the real number domain A region in , and let is an analytical function with the following integral form, the path And function With interpolation nodes The interpolation polynomial p of f m-1 for: Here w m ψ(z) = (z - θ1)…(z - θ m ); Next, we give the error integral expression of the Arnoldi approximation f(A)b: Since the matrix spec(A)∈Ω, and the vector If known, use f m represents the mth Arnoldi approximation of f(A)b; Because p m-1 The interpolation nodes are θ1, θ2, …, θ m The interpolation polynomial of f is obtained by subtracting p from expression (11). m-1 ,get: Substitute the matrix A for z in formula (10), and then multiply it by the vector b; since p m-1 (A) b = f m ,get: Here The above formula shows that the Arnoldi approximation f m The error with f(A)b is an error matrix function e m (A) Multiplication with a vector.

4. The transient electromagnetic fast forward modeling method based on the integral restart algorithm according to claim 3 is characterized in that: The step 4 is specifically as follows: For the integration node t i ∈Γ and weight factor ω i , e in formula (11) m Any numerical integration formula for (A) has the following form Let the integration nodes and weight factors in formula (12) be fixed in all restart cycles; in the first restart cycle, that is, w m (t i )=1, we can get the first Arnoldi approximation: In subsequent restart processes, i.e., k ≥ 2, the approximate value of each restart cycle can be calculated by using the error expression (12). Finally, the following expression is obtained 5. The transient electromagnetic fast forward modeling method based on the integral restart algorithm according to claim 4 is characterized in that: The step five is specifically as follows: The Arnoldi approximation for each restart cycle calculated in step 4 Thus, the corresponding matrix exponential function f(A)b is calculated. Since f(A)=exp(-tA), b=B0, and the Arnoldi algorithm is essentially equivalent to an iterative algorithm that adaptively selects the time step, in order to speed up the calculation, formula (7) is rewritten as: B(t i )=exp(-h i A)B i-1 (18) Here, t i =t i-1 +h i ,i=1,2,…,nt,H i is the time step, nt is the number of discrete time channels; Then, the geoelectric model of step 1 is forward simulated using the rewritten formula (18) and GPU parallel processing. In the GPU parallel processing, the specific processing steps are as follows: first, the matrix data is copied from the main memory to the GPU memory, the calculation process is guided by the CPU, and then intensive numerical operations are quickly performed on each core of the GPU. Finally, the calculation results are returned to the main memory to complete the forward modeling process.

Images