Octree grid-based semi-airborne full-time domain electromagnetic forward method

Abstract

The application discloses a semi-airborne full-time electromagnetic forward method based on an Octree grid, which performs spatial discretization on time-domain Maxwell equations, thereby effectively reducing unknowns of equations to be solved. For full-time electromagnetic field numerical simulation, the application combines the high calculation precision of the second-order backward Euler method with the high OFF-Time calculation efficiency of the model reduction algorithm, adopts the second-order BDF and the model reduction algorithm to respectively solve ON-Time and OFF-Time time-domain control equations, forms a joint algorithm to realize simulation of semi-airborne time-domain electromagnetic full-time response, and finally completes the forward calculation process. Through numerical simulation verification, in the full-time transient electromagnetic forward process, the application can keep small occupied memory and also has higher calculation efficiency, and the calculation precision is not lower than that of other methods.

Description

Technical Field

[0001] This invention belongs to the field of transient electromagnetic detection, specifically a semi-airborne full-time-domain electromagnetic forward modeling method based on the Octree grid. Background Technology

[0002] Time-domain airborne electromagnetic surveying is a major category of airborne electromagnetic surveying. It involves using an airborne electromagnetic instrument mounted on an aircraft (or helicopter sling) to transmit pulsed currents of a given interval and shape to a large transmitting coil mounted on the aircraft (or helicopter sling), emitting a strong primary electromagnetic field into the ground. A receiving coil mounted in the aircraft (or helicopter pod) measures the attenuation characteristics of the secondary electromagnetic field generated by the underground geological body in different directions and at different times during the transmission interruption (and during transmission). Transient electromagnetic forward modeling is then used to obtain the desired internal geological structure. Semi-airborne time-domain electromagnetic surveying is one type of time-domain airborne electromagnetic surveying. This method is significantly affected by the transmitted waveform; therefore, the forward modeling must consider the transmitted waveform. However, traditional full-time three-dimensional numerical simulation of time-domain electromagnetic surveying is computationally expensive (due to the large amount of computation). In other words, current processing methods require a large amount of memory for good computational efficiency, and very low computational efficiency for small memory usage.

[0003] Therefore, the research direction required by this invention is to provide a new method that can achieve high computational efficiency while maintaining a small memory footprint during semi-airborne electromagnetic full-time forward modeling. Summary of the Invention

[0004] To address the problems existing in the prior art, this invention provides a semi-aircraft full-time electromagnetic forward modeling method based on Octree grids, which can achieve high computational efficiency while maintaining a small memory footprint during full-time transient electromagnetic forward modeling.

[0005] To achieve the above objectives, the technical solution adopted by this invention is: a semi-aircraft full-time-domain electromagnetic forward modeling method based on Octree grids, the specific steps of which are as follows:

[0006] Step 1: Model parameterization: Based on known geological data, establish a geophysical geoelectric model;

[0007] Step 2, Spatial Discretization: The time-domain Maxwell equations are spatially discretized using the Octree mesh-based pseudo-finite volume method to obtain the spatially discretized matrix expressions of the governing equations.

[0008] Step 3: Obtain ON-Time electromagnetic field data: Based on the matrix expression obtained in Step 2 and the backward Euler method, solve for the ON-Time electromagnetic field data.

[0009] Step 4: Obtain OFF-Time Electromagnetic Field Data: Based on the matrix expression obtained in Step 2 and the model reduction method, the OFF-Time electromagnetic field data is obtained by solving the problem.

[0010] 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.

[0011] Furthermore, step two specifically involves:

[0012] The time-domain electromagnetic full-time response can be obtained by solving Maxwell's equations, neglecting displacement current, as follows:

[0013]

[0014] Where e and b are the electric field strength and magnetic induction intensity, respectively, t represents time, μ0 is the permeability of free space, σ 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, a first-kind boundary condition is adopted:

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

[0016] Using the pseudo-finite volume method to discretize the Maxwell equations in the time domain, by introducing test parameters w and f that lie in the same Sobolev space as e and b, and then productting them with equation (1), we obtain:

[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 the electromagnetic field on the Octree grid follows the Yee spatial discretization scheme. The electric field e is defined at the midpoint of the edge, and its three components are e0, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e2 ... x ,e y and e z The direction is parallel to the edge. Simultaneously, 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 of the normal vector is parallel to the mesh surface;

[0021] By discretizing the curl and the inner product, equation (3) can be written in the following form:

[0022]

[0023] Since w and f are arbitrary test parameters, w and f can be eliminated from both sides of equation (5), resulting in a discrete matrix expression for the control equation space:

[0024]

[0025] Furthermore, step three specifically includes:

[0026] After eliminating the electric field e from equation (6), for Discretization using the second-order backward Euler formula yields the governing equation for the magnetic flux density b:

[0027]

[0028] Here, the coefficient matrix is ​​denoted. t k =t k-1 +Δt k k = 1, 2, ..., n, Δt k It is the ON-Time time step and t n =t off At the initial time t0 = 0, the source term is 0, i.e., s0 = s -1 =0, therefore 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 ON-Time electromagnetic field data can be obtained.

[0029] Furthermore, step four specifically includes:

[0030] The governing equation for the magnetic flux density b after the current is turned off is as follows:

[0031]

[0032] Let matrix For solving 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, the direct model reduction is adopted for the governing equation, and the residual threshold tol is used to control the forward modeling process, so as to obtain the off-time electromagnetic field data.

[0033] Furthermore, the specific process of directly solving the OFF-Time electromagnetic field data using a reduced-order model is as follows:

[0034] First, construct the matrix I+γA off The projection subspace of vector u0:

[0035]

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

[0037] basis vector matrix V m sum matrix H m It has the following relationship:

[0038]

[0039] in It is a matrix composed of a series of orthogonal basis vectors, satisfying N is the number of unknowns to be solved. Formula (10) can be rewritten as follows:

[0040]

[0041] Here From formula (10):

[0042]

[0043] Using the Arnoldi matrix Multiplying both sides of equation (8) simultaneously yields:

[0044]

[0045] By reducing the order of the model of equation (8) and letting The following low-order partial differential equations are obtained:

[0046]

[0047] Observing the above equation, we find 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), we obtain the magnetic induction intensity and induced electromotive force of the original control equation:

[0048]

[0049] This allows us to obtain the electromagnetic field data for the off-time.

[0050] Compared with existing technologies, this invention employs a pseudo-finite volume method based on Octree grids to spatially discretize the Maxwell equations in the time domain. Compared with existing Tensor grids, Octree grids have significant advantages in local refinement, avoiding redundant grids at the boundaries due to local refinement, thus effectively reducing the unknowns in the equations to be solved. For full-time electromagnetic field numerical simulation, both the full-time back-Eulerian method and the model order reduction method have been proven effective. However, the former requires a large number of matrix decompositions and hundreds or thousands of matrix back-substitution processes, while the latter requires multiple constructions of the projected subspace during the ON-Time period, and its computational accuracy is lower than that of the second-order back-Eulerian method. To solve the above problems, this invention combines the high computational accuracy of the second-order back-Eulerian method with the high computational efficiency of the SAI-Krylov algorithm in the OFF-Time period. It uses the second-order BDF and the SAI-Krylov algorithm to solve the time domain governing equations in the ON-Time and OFF-Time periods, forming a joint algorithm to simulate the full-time electromagnetic response in the semi-aircraft time domain. Numerical simulations have verified that this invention can achieve high computational efficiency while maintaining a small memory footprint during the real-time transient electromagnetic forward modeling process; and its computational accuracy is basically the same as other existing methods that have a large memory footprint or long computation time. Attached Figure Description

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

[0052] In the figure: (a) magnetic field strength is assigned to the center of the grid cell; (b) electric field strength is assigned to the center of the edge.

[0053] Figure 2 It is a schematic diagram of any type of transmission waveform.

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

[0055] In the figure: (a) transmitted waveform; (b) time step.

[0056] Figure 4 These are the results of forward modeling the HELTITEM MULTIPULSE system using four different methods.

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

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

[0059] In the figure: (a) transmitted waveform; (b) time step.

[0060] Figure 6 These are the results of forward modeling the VTEM system using four different methods.

[0061] In the diagram: (a)b z (b)db z / dt. Detailed Implementation

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

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

[0064] Step 1: Model parameterization: Based on known geological data, establish a geophysical geoelectric model.

[0065] Step 2: Spatial Discretization: The time-domain Maxwell equations are spatially discretized using the Octree grid-based pseudo-finite volume method to obtain the spatially discretized matrix expressions of the governing equations, specifically:

[0066] The time-domain electromagnetic full-time response can be obtained by solving Maxwell's equations, neglecting displacement current, as follows:

[0067]

[0068] Where e and b are the electric field strength and magnetic induction intensity, respectively, t represents time, μ0 is the permeability of free space, σ 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, a first-kind boundary condition is adopted:

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

[0070] Using the pseudo-finite volume method to discretize the Maxwell equations in the time domain, by introducing test parameters w and f that lie in the same Sobolev space as e and b, and then productting them with equation (1), we obtain:

[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 an example, the inner product operation is defined as follows:

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

[0074] The space is discretized using an Octree grid based on orthogonal regular rectangles. The distribution of the electromagnetic field on the Octree grid follows the Yee spatial discretization scheme. The electric field e is defined at the midpoint of the edge, and its three components are e0, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e2 ...x ,e y and e z The direction is parallel to the edge. Simultaneously, 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 , which is a normal vector whose direction is parallel to the mesh surface.

[0075] By discretizing the curl and the inner product, equation (3) can be written in the following form:

[0076]

[0077] Since w and f are arbitrary test parameters, w and f can be eliminated from both sides of equation (5), resulting in a discrete matrix expression for the control equation space:

[0078]

[0079] Step 3: Obtain ON-Time Electromagnetic Field Data: Based on the matrix expression obtained in Step 2 and the backward Euler method, the ON-Time electromagnetic field data is obtained by solving the problem. Specifically:

[0080] After eliminating the electric field e from equation (6), for Discretization using the second-order backward Euler formula yields the governing equation for the magnetic flux density b:

[0081]

[0082] Here, the coefficient matrix is ​​denoted. t k =t k-1 +Δt k k = 1, 2, ..., n, Δt k It is the ON-Time time step and t n =t off ; Figure 2 This is a schematic diagram of an arbitrary transmitted waveform. At the initial time t0 = 0, the source term is 0, i.e., s0 = s -1 =0, therefore the electromagnetic field is also zero, e0 = e -1 =0, b0=b -1 =0; Observing equation (7), it can be seen that since the coefficient matrix changes with the time step, if we want to reduce the number of matrix decompositions, we can use linear equal interval discretization of time during the ON-Time period as much as possible; by using the direct solver PARDISO to iteratively solve equation (7), we can obtain the ON-Time electromagnetic field data.

[0083] Step 4: Obtain Off-Time Electromagnetic Field Data: Based on the matrix expression obtained in Step 2 and the model reduction method, the off-time electromagnetic field data is obtained by solving the problem. Specifically:

[0084] The governing equation for the magnetic flux density b after the current is turned off is as follows:

[0085]

[0086] Let matrix For solving equation (8), the model reduction method is a very efficient solution technique. However, unlike the model reduction of the analytical expression of the field, the direct model reduction is used for the governing equation, and the residual threshold is used to control the forward modeling process. The specific process is as follows:

[0087] First, construct the matrix I+γA off The projection subspace of vector u0:

[0088]

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

[0090] basis vector matrix V m sum matrix H m It has the following relationship:

[0091]

[0092] in It is a matrix composed of a series of orthogonal basis vectors, satisfying N is the number of unknowns to be solved. Formula (10) can be rewritten as follows:

[0093]

[0094] Here From formula (10):

[0095]

[0096] Using the Arnoldi matrix Multiplying both sides of equation (8) simultaneously yields:

[0097]

[0098] By reducing the order of the model of equation (8) and letting The following low-order partial differential equations are obtained:

[0099]

[0100] Observing the above equation, we find 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), we obtain the magnetic induction intensity and induced electromotive force of the original control equation:

[0101]

[0102] This allows us to obtain the electromagnetic field data for the off-time.

[0103] 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.

[0104] Numerical simulation verification:

[0105] Two typical emission waveforms were selected to excite the electromagnetic field. The forward modeling process was performed using the algorithm of this invention (abbreviated as BDF-MOROctree), as well as existing Tensor mesh BDF algorithms (abbreviated as BDF Tensor), BDF-MOR algorithms (abbreviated as BDF-MOR Tensor), and Octree algorithms (abbreviated as BDF Octree), as detailed below:

[0106] 1. Electromagnetic response simulation of the Helium Multiplyse system

[0107] First, the theoretical current waveform of the HELTITEM MULTIPULSE system is used, as shown below. Figure 3 As shown in Figure a, the waveform has a pulse width of 18ms, consisting of a 4ms half-sine wave followed by a 10.5ms off-time, then a 1ms trapezoidal wave followed by a 2.5ms off-time. The BDF algorithm uses linearly equal-interval discretization during the ON-Time periods of the half-sine and trapezoidal waves, with time steps of 40µs and 10µs respectively. During their corresponding OFF-Time periods, the time steps are set to 2... k-1 ×10 -7 s, k = 1, 2, ... gradually increase as follows Figure 3 As shown in b. The algorithm of this invention uses a model order reduction algorithm to quickly solve the problem within 10.5ms after the half-sine wave is turned off and within 2.5ms after the trapezoidal wave is turned off.

[0108] Table 1 shows the memory consumption and forward modeling time for four solution methods for the MULTIPULSE system. It can be seen that Octree mesh provides better mesh refinement compared to Tensor mesh, effectively avoiding the generation of a large number of redundant meshes. From a time cost perspective, Octree mesh generates fewer unknowns than Tensor mesh, resulting in lower time costs. Furthermore, the algorithm of this invention requires only one matrix factorization in the off-time using the model order reduction method, further improving the computational efficiency of the full-time electromagnetic field.

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

[0110]

[0111] Figure 4 a and Figure 4 b represents the MULTIPULSE system, which is always-on. z and db z / dt shows 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. Under this premise, it can be concluded that the present invention not only consumes less memory in these forward modeling methods, but also maintains a very short forward modeling time. Therefore, it proves that under the waveform electric field generated by the system, the method of the present invention has advantages over the existing methods in terms of memory consumption and forward modeling time.

[0112] 2. Electromagnetic response simulation of VTEM system

[0113] To further verify the ability of the algorithm of this invention to simulate the full-time electromagnetic response of complex waveforms, the actual transmitted waveform of the VTEM system was used, as shown in the figure. Figure 5 As shown in a. From Figure 5 As can be seen, the VTEM waveform is not smooth, exhibiting a sawtooth pattern during the ON-time period. To fully fit the current waveform, a linearly discrete time step of 40µs is used to discretize the VTEM waveform from 0ms to 8ms. For the electromagnetic response from 8ms to 16ms, the BDF algorithm still uses a varying time step for iterative solution. Figure 5 As shown in b, the present invention uses a model order reduction algorithm to solve the problem quickly.

[0114] Table 2 shows the memory consumption and forward modeling time for four solution methods for the VTEM system. Since the geological anomaly and mesh subdivision remain unchanged, the memory consumption is the same as the data in Table 1. Compared with the forward modeling time of the other three methods, the algorithm of this invention (i.e., BDF-MOR Octree) has a speedup of 46.5, 9.9, and 4.8, respectively.

[0115] Table 2 Memory consumption and forward modeling time of VETM transmitted waveforms

[0116]

[0117]

[0118] Figure 6 a and Figure 6 b represents the full-time b of the VTEM system. z and db z / dt, as can be seen from the figure, the b of the four methods z and db z / dt matches well across all time periods, db z / dt exhibits lower consistency only at the non-differentiable points of the current waveform. This indicates that the accuracy of the four methods is essentially the same under this waveform condition. Under this premise, it can be concluded that the present invention not only consumes less memory among these forward modeling methods, but also maintains a very short forward modeling time.

[0119] The experimental verification of two different waveforms shows that, under different waveform electric field conditions, the present invention, while maintaining the same accuracy as the existing method, not only consumes less memory but also has a shorter forward modeling time (i.e., higher computational efficiency) compared to the existing method.

[0120] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A semi-aircraft full-time electromagnetic forward modeling method based on Octree grids, characterized in that, The specific steps are as follows: Step 1: Model parameterization: Based on geological data, establish a geophysical geoelectric model; Step 2: Spatial Discretization: The time-domain Maxwell equations are spatially discretized using the Octree grid-based pseudo-finite volume method to obtain the spatially discretized matrix expressions of the governing equations. Step 3: Obtain ON-Time Electromagnetic Field Data: Based on the matrix expression obtained in Step 2 and the backward Euler method, the ON-Time electromagnetic field data is obtained by solving the problem. Specifically: Eliminate the electric field from equation (6) Afterwards, The magnetic flux density is obtained by discretization using the second-order backward Euler formula. The governing equations: Here, the coefficient matrix is ​​denoted. , , , It is the time step of ON-Time and At the initial moment When the source term is 0, that is Therefore, the electromagnetic field is also zero. , The electromagnetic field data in ON-Time were obtained by iteratively solving equation (7) using the direct solver PARDISO. Step 4: Obtain OFF-Time Electromagnetic Field Data: Based on the matrix expression obtained in Step 2 and the model reduction method, the OFF-Time electromagnetic field data is obtained by solving the problem. 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. The semi-aircraft full-time electromagnetic forward modeling method based on Octree grids according to claim 1, characterized in that, Step two specifically involves: The time-domain electromagnetic full-time response is obtained by solving Maxwell's equations, neglecting displacement current, as follows: in and These are electric field strength and magnetic induction intensity, Indicates time, It is the vacuum permeability. Indicates the conductivity of different media. This is the source term; to ensure that the electromagnetic field decays to zero at a sufficiently far boundary, a first-type boundary condition is used: The Maxwell equations in the discrete-time domain are adopted using the pseudo-finite volume method, by introducing... and Test parameters located in the same Sobolev space and And by taking the inner product of equation (1), we get: The above express and The inner product, express and The inner product of the magnetic induction intensity and test parameters For example, the inner product operation is defined as follows: The space is discretized using an Octree grid based on orthogonal regular rectangles. The distribution of the electromagnetic field on the Octree grid follows the Yee spatial discretization scheme. Defined at the midpoint of the edge, let its three components be respectively and The direction is parallel to the edge; at the same time, the magnetic induction intensity Defined at the center of the mesh surface, let its three components be respectively and The direction of the normal vector is parallel to the mesh surface; By discretizing the curl and the inner product, equation (3) can be written in the following form: because and Since the test parameters are arbitrary, both sides of equation (5) can be eliminated. and This yields the spatially discrete matrix expression of the governing equations.

3. The semi-aircraft full-time electromagnetic forward modeling method based on Octree grids according to claim 1, characterized in that, Step four specifically involves: Magnetic flux density after current is turned off The governing equations are expressed as follows: Let matrix The governing equations are reduced in order using a direct model, and the forward modeling process is controlled by a residual threshold, thereby obtaining off-time electromagnetic field data.

4. The semi-aircraft full-time electromagnetic forward modeling method based on Octree grids according to claim 3, characterized in that, The specific process of directly solving the OFF-Time electromagnetic field data by reducing the order of the model is as follows: First, construct the matrix. sum vector Projected subspace: Here m is the order of the SAI-Krylov subspace. It is a displacement parameter. It is the identity matrix; basis vector matrix sum matrix It has the following relationship: in It is a matrix composed of a series of orthogonal basis vectors, satisfying N is the number of unknowns to be solved. Formula (10) can be rewritten as follows: Here From formula (10), we get: Using the Arnoldi matrix Multiplying both sides of equation (8) simultaneously yields: By reducing the order of the model of equation (8) and letting The following low-order partial differential equations are obtained: Observing the above equation, we find that the order of equation (14) is m, which is much smaller than the order N of the original governing equation (8), that is... The magnetic flux density and induced electromotive force of the original governing equations are obtained using formula (15): This allows us to obtain the electromagnetic field data for the off-time.

Images