Semi-aviation full-time-domain electromagnetic forward modeling method based on Octree grid

Abstract

The invention discloses a semi-aviation full-time-domain electromagnetic forward modeling method based on an Octree grid, and the method carries out the spatial discretization of a time domain Maxwell equation, thereby effectively reducing the unknown number of an equation to be solved. For full-time electromagnetic field numerical simulation, the advantages of high calculation precision of a second-order backward Euler method and high calculation efficiency of a model order reduction algorithm in OFF-Time are combined, time domain control equations of ON-Time and OFF-Time are solved respectively by adopting the second-order BDF and the model order reduction algorithm, and a joint algorithm is formed to realize simulation of semi-aviation time domain electromagnetic full-time response. And finally completing a forward modeling calculation process. Through numerical simulation verification, in the full-time transient electromagnetic forward modeling process, the calculation efficiency can be higher on the premise of keeping small occupied memory; and the calculation precision is not lower than that of other existing methods.

Description

Technical Field

[0001] The present invention belongs to the field of transient electromagnetic detection, and particularly relates to a semi-aerial full-time domain electromagnetic forward modeling method based on Octree grids. Background Art

[0002] Time-domain airborne electromagnetic survey is a major category in airborne electromagnetic surveys. An airborne electromagnetic instrument installed on an aircraft (or suspended by a helicopter) sends pulsed currents with a given interval and shape to a large transmitting coil installed on the aircraft (or suspended by a helicopter), emitting a strong primary electromagnetic field into the ground. The receiving coil installed in the aircraft (or helicopter) pod measures the characteristics of the secondary electromagnetic field generated by underground geological bodies decaying in different directions and at different times during the period when the transmission is disconnected (and during the transmission). Transient electromagnetic forward modeling is used to obtain the internal structure of the geological area to be detected. The semi-aerial time-domain electromagnetic field is one type of time-domain airborne electromagnetic survey. This method is significantly affected by the transmitted waveform, so the forward simulation must consider the transmitted waveform. However, the traditional full-time three-dimensional numerical simulation of time-domain electromagnetics has an excessively high computational cost (the reason is that a large amount of calculation is generated). That is, for the current processing method, if better computational efficiency is required, a large amount of memory needs to be occupied; if only a small amount of memory is occupied, the computational efficiency is very low.

[0003] Therefore, how to provide a new method that can maintain a small memory footprint and high computational efficiency during the semi-aerial electromagnetic full-time 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 semi-aerial full-time domain electromagnetic forward modeling method based on Octree grids, which can maintain a small memory footprint and high computational efficiency during the full-time transient electromagnetic forward modeling process.

[0005] To achieve the above object, the technical solution adopted by the present invention is as follows: A semi-aerial full-time domain electromagnetic forward modeling method based on Octree grids, 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 mimetic finite volume method based on Octree grids to perform spatial discretization on the time-domain Maxwell equation to obtain a matrix expression for the spatial discretization of the control equation;

[0008] Step 3. Obtain the electromagnetic field data of ON-Time: According to the matrix expression obtained in Step 2 and combined with the backward Euler method, the electromagnetic field data of ON-Time is solved and obtained;

[0009] Step 4. Obtain the electromagnetic field data of OFF-Time: According to the matrix expression obtained in Step 2 and combined with the model reduction method, the electromagnetic field data of OFF-Time is solved and obtained;

[0010] Step 5. Data integration: Integrate the electromagnetic field data of ON-Time obtained in Step 3 and the electromagnetic field data of OFF-Time obtained in Step 4, so as to complete the forward simulation of the detection area.

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

[0012] The full-time response of time-domain electromagnetic fields can be obtained by solving the Maxwell equations that neglect displacement currents. The equations are as follows:

[0013]

[0014] where e and b are the electric field strength and magnetic induction intensity respectively, t represents time, μ 0 is the vacuum magnetic permeability, σ represents the conductivity of different media, s is the source term; to ensure that the electromagnetic field decays to zero at a sufficiently far boundary, the first-kind boundary condition is adopted:

[0015] b×n = 0, n×e = 0 (2)

[0016] The mimetic finite volume method is used to discretize the time-domain Maxwell equations. By introducing the test parameters w and f that are in the same Sobolev space as e and b, and taking the inner product with Equation (1), we get:

[0017]

[0018] The above (b, f) represents the inner product of b and f, and (e, f) represents the inner product of e and f. Taking the magnetic induction intensity b and the test parameter f as an example, the inner product operation is defined as follows:

[0019] (b, f) = ∫ V b·fdV (4)

[0020] The space is discretized using an Octree grid based on orthogonal regular rectangles. The distribution of electromagnetic fields on the Octree grid follows the Yee space discretization scheme. The electric field e is defined at the midpoint of the edge, and its three components are respectively e x , e y and e z , and the direction is parallel to the edge. At the same time, the magnetic induction intensity b is defined at the center of the grid surface, and its three components are respectively bx , b y and b z , the direction is parallel to the normal vector of the grid plane;

[0021] Through curl discretization and inner product discretization, Equation (3) is written in the following form:

[0022]

[0023] Since w and f are arbitrarily introduced test parameters, w and f can be cancelled on both sides of Equation (5), obtaining the matrix expression for the spatial discretization of the governing equation:

[0024]

[0025] Furthermore, the specific content of Step 3 is as follows:

[0026] After eliminating the electric field e from Equation (6), for Adopting the second-order backward Euler formula for discretization, the governing equation for the magnetic induction intensity b is obtained:

[0027]

[0028] Here, the coefficient matrix is denoted as t k = t k-1 + Δt k , k = 1, 2, …, n, Δt k is the time step of ON-Time and t n = t off ; At the initial time t 0 = 0, the source term is 0, that is, s 0 = s -1 = 0, so the electromagnetic field is also zero, e 0 = e -1 = 0, b 0 = b -1 = 0; By using the direct solver PARDISO to iteratively solve Equation (7), the electromagnetic field data of ON-Time is obtained.

[0029] Furthermore, the specific content of Step 4 is as follows:

[0030] The governing equation for the magnetic induction intensity b after turning off the current is expressed as follows:

[0031]

[0032] Denote the matrix For the solution of equation (8), the model reduction method is a very efficient solution technique. However, it is different from the model reduction of the analytical expression of the field. Therefore, direct model reduction is applied to the governing equations, and the forward process is controlled by the residual threshold tol to obtain the electromagnetic field data of the OFF-Time.

[0033] Furthermore, the specific process of directly solving the electromagnetic field data of the OFF-Time by model reduction is as follows:

[0034] First, construct the projection subspace of the matrix I + γA off and the vector u 0 :

[0035]

[0036] Here, u 0 = A off b off , m is the order of the SAI-Krylov subspace, γ is the shift parameter, and I is the identity matrix;

[0037] The basis vector matrix V m and the matrix H m have the following relationship:

[0038]

[0039] where is a matrix composed of a series of orthogonal basis vectors, satisfying N is the number of unknowns to be solved, Rewrite formula (10) as the following formula:

[0040]

[0041] Here Obtained from formula (10):

[0042]

[0043] Multiply both sides of equation (8) by the Arnoldi matrix simultaneously to obtain:

[0044]

[0045] By reducing the order of equation (8) and setting the following low-order partial differential equation is obtained:

[0046]

[0047] Observing the above equation, it is found that the order of equation (14) is m, which is much smaller than the order N of the original control equation (8), i.e., m << N; through formula (15), the magnetic induction intensity and induced electromotive force of the original control equation are obtained:

[0048]

[0049] Thus, the electromagnetic field data of OFF-Time is solved and obtained.

[0050] Compared with the prior art, the present invention adopts the mimetic finite volume method based on Octree grid to perform spatial discretization on the time-domain Maxwell equation; compared with the existing Tensor grid, the Octree grid has significant advantages in local encryption, which can avoid generating redundant grids at the boundary due to local encryption, thereby effectively reducing the unknowns of the equation to be solved. For the full-time electromagnetic field numerical simulation, both the full-time backward Euler method and the model reduction method have been proven effective. However, the former requires a large number of matrix decompositions and hundreds of matrix back-substitution processes, and the latter needs to construct the projection subspace multiple times during ON-Time, and its calculation accuracy is lower than that of the second-order backward Euler method. To solve the above problems, the present invention combines the advantages of the high calculation accuracy of the second-order backward Euler method and the high calculation efficiency of the SAI-Krylov algorithm in OFF-Time, and adopts the second-order BDF and the SAI-Krylov algorithm to solve the time-domain control equations of ON-Time and OFF-Time respectively, forming a combined algorithm to realize the simulation of the full-time response of the semi-aerial time-domain electromagnetic field. Through numerical simulation verification, during the full-time transient electromagnetic forward modeling of the present invention, it can maintain a small memory occupancy and also have a relatively high calculation efficiency; and its calculation accuracy is basically the same as that of other existing methods with a large memory occupancy or a long calculation time. Brief Description of the Drawings

[0051] Figure 1 is a schematic diagram of the electromagnetic field distribution after discretization of the present invention.

[0052] In the figure: (a) The magnetic induction intensity is assigned to the center of the grid surface element; (b) The electric field intensity is assigned to the center of the edge.

[0053] Figure 2 is a schematic diagram of any emission waveform.

[0054] Figure 3 is a schematic diagram of the emission current parameters of the HELTITEM MULTIPULSE system in the numerical simulation verification.

[0055] In the figure: (a) Emission waveform; (b) Time step.

[0056] Figure 4These are the forward modeling result diagrams of the HELTITEM MULTIPULSE system using four methods respectively.

[0057] In the figure: (a) b z ; (b) db z / dt.

[0058] Figure 5 This is the schematic diagram of the transmitting current parameters of the VTEM system in the numerical simulation verification.

[0059] In the figure: (a) Transmitting waveform; (b) Time step.

[0060] Figure 6 These are the forward modeling result diagrams of the VTEM system using four methods respectively.

[0061] In the figure: (a) b z ; (b) db z / dt. Detailed implementation method

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

[0063] Example: The specific steps are as follows:

[0064] Step 1: Model parameterization: Establish a geophysical electro - geological model according to the known geological data.

[0065] Step 2: Spatial discretization processing: Use the mimetic finite - volume method based on the Octree grid to perform spatial discretization on the time - domain Maxwell equation to obtain the matrix expression of the spatial discretization of the control equation. Specifically:

[0066] The full - time response of the time - domain electromagnetic field can be obtained by solving the Maxwell equations neglecting the displacement current. The equations are as follows:

[0067]

[0068] Where e and b are the electric - field strength and magnetic - induction intensity respectively, t represents time, μ 0 is the vacuum permeability, σ represents the conductivity of different media, s is the source term; to ensure that the electromagnetic field decays to zero at a sufficiently far boundary, the first - type boundary condition is adopted:

[0069] b×n = 0, n×e = 0 (2)

[0070] Discretize the time - domain Maxwell equations using the mimetic finite - volume method. By introducing the test parameters w and f in the same Sobolev space as e and b and taking the inner product with equation (1), we get:

[0071]

[0072] The above (b, f) represents the inner product of b and f, and (e, f) represents the inner product of e and f. Taking the magnetic induction intensity b and the test parameter f as examples, the inner product operation is defined as follows:

[0073] (b, f) = ∫ V b·fdV (4)

[0074] Using an Octree grid based on orthogonal regular rectangles to discretize the space, the distribution of the electromagnetic field on the Octree grid follows the Yee space discretization scheme. The electric field e is defined at the midpoint of the edge, and its three components are respectively e x , e y and e z , with the direction parallel to the edge. At the same time, the magnetic induction intensity b is defined at the center of the grid surface, and its three components are respectively b x , b y and b z , with the direction parallel to the normal vector of the grid surface.

[0075] Through curl discretization and inner product discretization, Equation (3) is written in the following form:

[0076]

[0077] Since w and f are arbitrarily introduced test parameters, w and f can be cancelled on both sides of Equation (5), obtaining the matrix expression for the spatial discretization of the control equation:

[0078]

[0079] Step 3: Obtain the electromagnetic field data of ON-Time: According to the matrix expression obtained in Step 2 and combined with the backward Euler method, the electromagnetic field data of ON-Time is solved, specifically:

[0080] After eliminating the electric field e from Equation (6), the second-order backward Euler formula is used to discretize , obtaining the control equation for the magnetic induction intensity b:

[0081]

[0082] Here, the coefficient matrix is denoted as t k = t k-1 + Δt k , k = 1, 2,..., n, Δt k is the time step of ON-Time and t n = t off ; Figure 2 is a schematic diagram of an arbitrary emission waveform. At the initial time t 0 = 0, the source term is 0, that is, s0 = s -1 = 0, so the electromagnetic field is also zero, e 0 = e -1 = 0, b 0 = b -1 = 0; Observing Equation (7), it can be seen that since the coefficient matrix changes with the time step, if you want to reduce the number of matrix decompositions, you can discretize the time linearly at equal intervals as much as possible during the ON-Time; By using the direct solver PARDISO to iteratively solve Equation (7), the electromagnetic field data for the ON-Time can be obtained.

[0083] Step 4. Obtain the electromagnetic field data for the OFF-Time: According to the matrix expression obtained in Step 2 and combined with the model reduction method, the electromagnetic field data for the OFF-Time can be solved, specifically:

[0084] The control equation for the magnetic induction intensity b after turning off the current is expressed as follows:

[0085]

[0086] Denote the matrix For the solution of Equation (8), the model reduction method is a very efficient solution technique, but it is different from the model reduction of the analytical expression of the field. Therefore, direct model reduction is used for the control equation, and the forward process is controlled by the residual threshold. The specific process is as follows:

[0087] First, construct the projection subspace of the matrix I + γA off and the vector u 0 :

[0088]

[0089] Here u 0 = A off b off , m is the order of the SAI-Krylov subspace, γ is the displacement parameter, and I is the identity matrix;

[0090] The basis vector matrix V m and the matrix H m have the following relationship:

[0091]

[0092] where is a matrix composed of a series of orthogonal basis vectors, satisfying N is the number of unknowns to be solved, Rewrite Equation (10) as the following formula:

[0093]

[0094] Here obtained from Equation (10):

[0095]

[0096] Using the Arnoldi matrix Multiply both sides of Equation (8) simultaneously to obtain:

[0097]

[0098] By reducing the order of the model of Equation (8) and setting the following low-order partial differential equation is obtained:

[0099]

[0100] Observing the above equation, it is found that the order of Equation (14) is m, which is much smaller than the order N of the original control equation (8), that is, m << N; through Equation (15), the magnetic induction intensity and induced electromotive force of the original control equation are obtained:

[0101]

[0102] Thus, the electromagnetic field data of OFF-Time is solved and obtained.

[0103] Step Five: Data integration: Integrate the electromagnetic field data of ON-Time obtained in Step Three and the electromagnetic field data of OFF-Time obtained in Step Four, so as to complete the forward simulation of the detection area.

[0104] Numerical simulation verification:

[0105] Select two typical emission waveforms to excite the electromagnetic field, and use the algorithm of the present invention (abbreviated as BDF-MOROctree), as well as the existing BDF algorithm of Tensor mesh (abbreviated as BDF Tensor), BDF-MOR algorithm (abbreviated as BDF-MOR Tensor) and BDF algorithm of Octree mesh (abbreviated as BDF Octree), respectively perform the forward process, specifically as follows:

[0106] 1. Electromagnetic response simulation of HELTITEM MULTIPULSE system

[0107] First, use the theoretical current waveform of the HELTITEM MULTIPULSE system as Figure 3As shown in a, the pulse width of this waveform is 18 ms, which successively includes a half-sine wave with a pulse width of 4 ms and an off period of 10.5 ms, followed by a trapezoidal wave with a pulse width of 1 ms and an off period of 2.5 ms. The BDF algorithm adopts linear equally spaced discretization during the ON-Time of the half-sine wave and the trapezoidal wave, and the time steps are 40 us and 10 us respectively. During the corresponding OFF-Time, the time step is multiplied by 2 k-1 ×10 -7 s, k = 1, 2, … and gradually increases as shown in Figure 3 b. The algorithm of the present invention uses the model reduction algorithm to quickly solve within 10.5 ms after the half-sine wave is turned off and within 2.5 ms after the trapezoidal wave is turned off.

[0108] Table 1 shows the memory consumption and forward modeling time of four solution methods for the MULTIPULSE system. It can be seen that the Octree mesh has a better encryption effect compared to the Tensor mesh and can effectively avoid generating a large number of redundant meshes. In terms of time cost, since the Octree mesh generates fewer unknowns compared to the Tensor mesh, the corresponding time cost is lower. In addition, the algorithm of the present invention only requires 1 matrix decomposition using the model reduction method during the OFF-Time, further improving the calculation efficiency of the full-time electromagnetic field.

[0109] Table 1 Memory consumption and forward modeling time of the MULTIPULSE waveform

[0110]

[0111] Figure 4 a and Figure 4 b are respectively the b z and db z / dt of the full-time of the MULTIPULSE system. It can be seen that the calculation results of the four methods are in good agreement, indicating that the forward modeling accuracy of each method is basically the same. On this premise, it can be concluded that the present invention not only has less memory consumption but also can maintain a very short forward modeling time among these forward modeling methods. Therefore, it is proved that the method of the present invention has advantages in terms of memory consumption and forward modeling time compared to the existing methods in the case of the waveform electric field generated by this system.

[0112] 2. Electromagnetic response simulation of the VTEM system

[0113] In order to further verify the ability of the algorithm of the present invention to simulate the full-time electromagnetic response of complex waveforms, the actual emission waveform of the VTEM system is used as shown in Figure 5 a. From Figure 5It can be seen that the VTEM waveform is not smooth and is serrated during the ON-Time. To fully fit the current waveform, the VTEM waveform from 0 ms to 8 ms is linearly discretized at equal intervals with a time step of 40 us. For the electromagnetic response from 8 ms to 16 ms, the BDF algorithm still uses a variable time step for iterative solution as Figure 5 shown in b, while the present invention uses a model reduction algorithm for rapid solution.

[0114] Table 2 shows the memory consumption and forward modeling time for four solution methods of the VTEM system. Since the geological anomalies and grid meshing remain unchanged, the memory consumption is the same as the data in Table 1. Comparing the forward modeling time of the algorithm of the present invention (i.e., BDF-MOR Octree) with that of the other three methods, the speedup ratios are 46.5, 9.9, and 4.8 respectively.

[0115] Table 2 Memory consumption and forward modeling time of the VETM emission waveform

[0116]

[0117]

[0118] Figure 6 a and Figure 6 b are the b of the full time of the VTEM system respectively z and db z / dt. It can be seen from the figure that the b z and db z / dt of the four methods are in good agreement throughout the full time period, and the consistency of db z / dt is only low at the non-differentiable points of the current waveform. This shows that the accuracies of the four methods are basically the same in this waveform case. On this premise, it can be concluded that the present invention not only has less memory consumption but also can maintain a very short forward modeling time among these forward modeling methods.

[0119] Based on the experimental verification of the two different waveforms, it can be seen that in the case of different waveform electric fields, on the premise of maintaining basically the same accuracy as the existing methods, the present invention not only has less memory consumption but also has a shorter forward modeling time (i.e., higher computational efficiency) compared with the existing methods.

[0120] The above are only the preferred embodiments of the present invention. It should be noted 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 semi-aeronautical full-time domain electromagnetic forward modeling method based on Octree grid, 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 pseudo-finite volume method based on the Octree grid to spatially discretize the Maxwell equations in the time domain and obtain the matrix expression of the spatial discretization of the control equations; Step 3, obtaining the electromagnetic field data of ON-Time: combining the matrix expression obtained in step 2 with the backward Euler method, thereby solving and obtaining the electromagnetic field data of ON-Time; Step 4, obtaining the electromagnetic field data of OFF-Time: According to the matrix expression obtained in step 2, combined with the model reduction method, the electromagnetic field data of OFF-Time is obtained; Step 5: Data integration: Integrate the ON-Time electromagnetic field data obtained in step 3 and the OFF-Time electromagnetic field data obtained in step 4 to complete the forward modeling of the detection area.

2. According to the Octree grid-based semi-aerial electromagnetic full-time forward modeling method of claim 1, it is characterized in that: The step 2 is specifically as follows: The full-time electromagnetic response in the time domain is obtained by solving the Maxwell equations ignoring the displacement current. The equations are as follows: Where e and b are the electric field intensity and magnetic induction intensity respectively, t represents time, μ0 is the vacuum magnetic permeability, σ represents the conductivity of different media, and s is the source term; to ensure that the electromagnetic field decays to zero at a sufficiently far boundary, the first type of boundary condition is adopted: b×n=0,n×e=0 (2) The pseudo-finite volume method is used to discretize the Maxwell equations in the time domain. By introducing the test parameters w and f in the same Sobolev space as e and b, and inner product with equation (1), we can obtain: The above (b, f) represents the inner product of b and f, and (e, f) represents the inner product of e and f. Taking the magnetic induction intensity b and the test parameter f as an example, the inner product operation is defined as follows: (b, f)=∫ V b·fdV (4) The Octree grid based on orthogonal regular rectangles is used to discretize the space. The distribution of the electromagnetic field on the Octree grid follows the Yee space discretization scheme. The electric field e is defined at the midpoint of the edge, and its three components are e x , e y and e z , the direction is parallel to the edge; at the same time, the magnetic induction intensity b is defined at the center of the grid surface, and its three components are b x , b y and b z , the direction is parallel to the normal vector of the mesh surface; Through the discretization of curl and inner product, equation (3) can be written as follows: Since w and f are arbitrary test parameters introduced, w and f can be eliminated on both sides of equation (5), and the matrix expression of the spatial discretization of the control equation is obtained:

3. The semi-aeroelectromagnetic full-time forward modeling method based on Octree grid according to claim 2 is characterized in that: The step three is specifically as follows: After eliminating the electric field e from equation (6), The second-order backward Euler formula is used to discretize the magnetic induction intensity b, and the control equation is obtained: Here we write the coefficient matrix t k =t k-1 +Δt k , k = 1, 2, ..., n, Δt k is the time step of ON-Time and t n =t off ; At the initial time t0 = 0, the source term is 0, that is, s0 = s -1 =0, so the electromagnetic field is also zero, e0=e-1=0, b0=b-1=0; by using the direct solver PARDISO to iteratively solve equation (7), the electromagnetic field data of ON-Time is obtained.

4. The semi-aeroelectromagnetic full-time forward modeling method based on Octree grid according to claim 3 is characterized in that: The step 4 is specifically as follows: The control equation of the magnetic induction intensity b after the current is turned off is expressed as follows: Remember the matrix Direct model reduction is used for the control equation, and the residual threshold is used to control the forward modeling process to obtain the OFF-Time electromagnetic field data.

5. The semi-aeroelectromagnetic full-time forward modeling method based on Octree grid according to claim 4 is characterized in that: The specific process of directly reducing the model to solve the electromagnetic field data of OFF-Time is as follows: First, construct the matrix I+γA off And the projection subspace of the vector u0: Here u0=A off b off , m is the order of SAI-Krylov subspace, γ is the shift parameter, and I is the identity matrix; Basis vector matrix V m and the matrix H m Has the following relationship: in is a matrix consisting of a series of orthogonal basis vectors that satisfies N is the number of unknowns to be solved, Rewrite formula (10) as follows: Here From formula (10), we get: Using the Arnoldi Matrix Multiplying both ends of equation (8) yields: By reducing the order of the model in equation (8), and setting The following low-order partial differential equation is obtained: By observing the above formula, it is found that the order of equation (14) is m, which is much smaller than the order N of the original control equation (8), that is, m<<N. By formula (15), the magnetic induction intensity and induced electromotive force of the original control equation are obtained: Thus, the electromagnetic field data of OFF-Time is obtained.

Images