Large loop source transient electromagnetic fast forward modeling method based on GPU acceleration

Abstract

The invention discloses a large-loop source transient electromagnetic fast forward modeling method based on GPU acceleration, and the method comprises the steps: firstly carrying out the parameter setting and preprocessing, then segmenting a large-loop source into a north edge, a south edge, an east edge and a west edge, and carrying out the discretization of each edge into N electric dipoles through employing a Gaussian-Legendre quadrature method; constructing a loop edge-frequency sampling point-wave number-time channel four-level parallel computing architecture to compute each edge; wherein each side in the GPU adopts a frequency-wave number-time parallel computing architecture, and parallel computing is performed on each electric dipole, so that the frequency-time domain response of each electric dipole is obtained; and finally, after calculation of all the electric dipoles is completed, vectorization calculation and summation are carried out on transient responses of all the electric dipoles on the GPU, a large loop source total field response at the current observation point is obtained, and a final time domain electromagnetic field value is output. According to the method, the advantages of large-scale parallel computing of the GPU can be efficiently and fully played, so that the computing efficiency is effectively improved on the premise of ensuring the computing precision.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of detection inversion calculation, and particularly relates to a GPU accelerated large loop source transient electromagnetic fast forward method. BACKGROUND

[0002] The transient electromagnetic method (TEM) is an important geophysical exploration method, which is widely used in mineral exploration, hydrogeology, engineering and environmental exploration and other fields. Among them, the large loop source device is frequently used due to its large detection depth and strong signal. The core of the forward simulation is to calculate the electromagnetic response of the underground medium under the excitation of the transmitted current. Although the one-dimensional forward model can provide an accurate solution in theory, the numerical calculation efficiency is a significant bottleneck in practical application. Especially in the inversion interpretation, a large number of forward calculations are needed, and the calculation speed directly determines the overall efficiency.

[0003] At present, the traditional large loop source transient electromagnetic one-dimensional forward simulation is usually executed in series on a central processing unit (CPU), and the main technical route is as follows: first, the large loop source is discretized into N electric dipole units, the Hankel transform is applied to calculate the electromagnetic field response of each dipole in the frequency domain; then the forward and cosine numerical filtering algorithm is used to convert the frequency domain response to the time domain, and finally the discrete electric dipole response is superimposed by the numerical integration method to obtain the time domain electromagnetic field response of the large loop source. However, the above traditional method has the following technical bottlenecks: 1. Low calculation efficiency: The large loop source forward usually divides the loop into multiple electric dipoles, and uses the Gauss-Legendre quadrature method to calculate the contribution of each segment and then superimposes. The electric dipole response calculation involves Hankel transform and sine-cosine transform, and 140-point Hankel transform coefficients and 250-point sine-cosine filter coefficients are commonly used for frequency-time domain conversion, which has a huge amount of calculation. Although the traditional CPU serial calculation method can realize the forward, the calculation speed is slow, especially when the time channel increases or the model is complex, it is difficult to meet the demand for large-scale data simulation in actual production or fast forward in actual exploration.

[0004] 2. Low memory usage efficiency: The repeated data transmission and storage operation in the forward calculation process reduces the overall performance.

[0005] 3. Conflict between precision and speed: In the existing technology, increasing the number of filter coefficients can improve the calculation precision, but it will significantly reduce the calculation speed; while reducing the number of filter coefficients can improve the efficiency, but it will lead to precision loss, especially in the calculation of late time channels.

[0006] In the prior art, although there is a method of synthesizing a loop source by a dipole (such as integrating an electric dipole along a loop), the calculation process is not fully optimized, and the parallel computing capability of a GPU is not fully utilized, resulting in insufficient efficiency in large-scale forward simulation or inversion solving. Therefore, in view of the complex recursive calculation, multiple integration and transformation and other mathematical processes in transient electromagnetic forward calculation, how to provide a new fast forward framework to build a specific algorithm architecture that can efficiently and fully utilize the large-scale parallel computing advantage of a GPU, so as to effectively improve the calculation efficiency while ensuring the calculation accuracy, is the research direction of the present application. SUMMARY

[0007] In view of the problems in the prior art, the present application provides a GPU accelerated large loop source transient electromagnetic fast forward method, which builds a multi-level parallel computing architecture that can efficiently and fully utilize the large-scale parallel computing advantage of a GPU, thereby effectively improving the calculation efficiency while ensuring the calculation accuracy.

[0008] To achieve the above purpose, the technical scheme adopted by the present application is as follows: a GPU accelerated large loop source transient electromagnetic fast forward method, comprising the following steps: Step one, parameter setting and preprocessing: first, determine the number of required setting parameters, and assign values to each parameter according to the actual situation of the current observation point, if a parameter is missing in the actual situation, a default value is automatically assigned; then, use the persistent variable technology to pre-load the Hankel transform coefficients and the sine and cosine filter coefficients, avoid repeated loading and transmission, reduce the file I / O overhead, and transmit the above parameters and coefficients to the GPU memory, thereby reducing the CPU-GPU communication transmission; Step two, large loop source discretization and parallel task planning: the large loop source is divided into north, south, east and west four sides, each side is discretized into N electric dipoles (integration points) by using the Gauss-Legendre quadrature method; and a four-level parallel computing architecture is used to calculate the four sides, fully tapping the multi-core computing potential of the GPU; Step three, calculate the frequency-time domain response of each electric dipole: use the frequency-wave number-time parallel computing architecture in the GPU to perform parallel calculation on each electric dipole in step two, thereby obtaining the frequency-time domain response of each electric dipole; Step four, synthesizing the large loop source response and outputting: after the calculation of all electric dipoles in step three is completed, the transient responses of all electric dipoles are calculated and summed (Gaussian integration) on the GPU to obtain the total field response of the large loop source at the current observation point.

[0009] Further, the step one sets parameters including parameters of the layered earth model (resistivity ρ and thickness h of each layer), loop size (Lud is the length of up and down, Llr is the width of left and right), loop center coordinates (X0, Y0), measurement point coordinates (x, y), transmitting current I, and observation time channel t.

[0010] Further, the step two plans parallel tasks, specifically: A. Parallel calculation of loop edges: the four edge calculation tasks of the large loop source are distributed to different processing units (streaming multiprocessors) of the GPU, realizing parallel calculation of each edge at the macro task level; the calculation of each loop edge is completely independent, and high load balancing is realized through parallel loops.

[0011] B. Parallel layer of Gauss integral points: in the calculation of each loop edge, the Gauss-Legendre integration method is used to discretize each edge into N electric dipoles, and this calculation is performed in parallel at the GPU thread block level. This approach can fully utilize the advantages of the SIMT architecture of the GPU.

[0012] C. Parallel layer of Hankel transform: in the calculation of the frequency domain response of a single electric dipole, the calculation of 47 Hankel transform sampling points is distributed to GPU thread block level parallel processing; through an optimized memory access mode, the efficiency of high-frequency data access is ensured.

[0013] D. Parallel calculation layer of time channels: for frequency-time domain conversion calculation, the calculation tasks of different time channels are parallelized to realize fine-grained data level parallel processing; through parallel processing of time channels, the overall calculation time is significantly reduced.

[0014] Further, the step three is specifically: I. Frequency domain and wave number domain sampling: log-uniform sampling is performed in the complex frequency plane and wave number domain, respectively, to obtain multiple frequency sampling points ωi and wave number points λj.

[0015] II. Construction of three-dimensional parameter array: the parameters of the layered earth model, frequency, and wave number are expanded into a three-dimensional parameter array (number of layers × number of wave points × number of frequency points), and transmitted to the GPU; the purpose of this step is to reconstruct the entire calculation task into a highly regular data parallel problem.

[0016] III. Reconstruction of wave impedance recursive calculation: the traditional layer-by-layer recursive calculation is converted into matrix multiplication operation, and the frequency dispersion coefficient F(λ,ω) of the ground wave impedance is calculated combined with the three-dimensional parameter array; and the matrix calculation unit of the GPU is used to speed up the operation, that is, each thread block or thread of the GPU independently calculates the complete recursive process under a specific (λ,ω) combination, thereby converting the serial recursion into independent parallel recursion tasks; IV. Optimization of Hankel transform: using the dispersion coefficient F(λ, ω) calculated in step III, the Hankel transform is performed by a linear filtering algorithm and vectorized calculation to obtain the frequency domain electromagnetic field response H z (ω). V. Improved sine-cosine transform: using a non-uniform sampling strategy, high-density sampling is maintained in the key frequency band, and the traditional 250-point sine and cosine filter coefficients are optimized to 87 points. Combined with the cubic spline interpolation algorithm, the precision loss caused by the reduction of coefficients is compensated, so that the frequency domain response H z (ω) is converted to time domain response h z (t) by sine-cosine transform. Batch processing interpolation calculation is performed through this process to reduce GPU-CPU data transmission.

[0017] Further, the step IV adopts 47-point Hankel transform coefficients in the process of Hankel transform, which replaces the traditional 140-point Hankel transform coefficients, thereby reducing the amount of calculation while ensuring accuracy.

[0018] Further, after obtaining the total field response of the large loop source at the current detection point in step four, the response is transmitted from the GPU memory back to the host memory, and then the final time domain electromagnetic field value is output.

[0019] Compared with the prior art, the present application has the following advantages: 1. Extremely high calculation speedup: The present application adopts a four-layer parallel architecture of "loop side-frequency sampling point-wave number-time channel" which matches the streaming multi-processor of GPU, fully taps the multi-core computing potential of GPU, and achieves a performance improvement of several orders of magnitude.

[0020] 2. Good calculation accuracy: The present application dynamically selects the optimal filter coefficient configuration and compensates for the precision loss caused by the reduction of filter coefficients through numerical interpolation method, ensuring that the calculation efficiency is improved while the accuracy meets the engineering requirements, ensuring the numerical accuracy of the forward simulation and being suitable for complex geological conditions.

[0021] 3. Strong practicability: The efficient computing power of the present application makes the method widely applicable to transient electromagnetic (TEM) data fast processing, large-scale data inversion, high-density sampling simulation and other practical production scenarios, greatly improving the work efficiency.

[0022] 4. Significant cost-effectiveness: The hardware of the present application only needs existing GPU hardware to achieve high-performance computing, without additional customization, reducing the investment cost of high-performance computing devices. BRIEF DESCRIPTION OF DRAWINGS

[0023] Figure 1 is the overall flowchart of the embodiment of the present application.

[0024] Figure 2 is a schematic diagram of discretization of a large loop source into electric dipoles in an embodiment of the present application.

[0025] Figure 3 is a schematic diagram of a frequency-wavenumber-time parallel computing architecture adopted by each electric dipole in an embodiment of the present application.

[0026] Figure 4 is a schematic diagram of GPU thread grid planning in an embodiment of the present application.

[0027] Figure 5 is a comparison diagram of filter coefficients using cosine transform in an embodiment of the present application and a conventional method.

[0028] Figure 6 is a comparison and relative error diagram of uniform half-space numerical solution and analytical solution obtained by an embodiment of the present application and a conventional method.

[0029] Figure 7 is a comparison diagram of calculation results under different layered earth model parameters in an embodiment of the present application. DETAILED DESCRIPTION

[0030] The present application will be further described below.

[0031] The hardware environment adopted in the present embodiment includes a workstation or server equipped with an NVIDIA GPU, a video memory capacity ≥ 8 GB, and a system memory ≥ 8 GB. The software environment includes a CentOS 8 operating system (Linux-x86_64), installation of MATLAB 2025b, CUDA Toolkit 12.2 and above versions, and configuration of MATLAB GPU computing environment.

[0032] As shown in Figures 1 to 3 , the present embodiment includes the following steps: Step one, parameter setting and pretreatment: first determine the number of required setting parameters, and assign values to each parameter according to the actual situation of the current observation point, including the parameters of the layered earth model (resistivity ρ and thickness h of each layer), loop size (Lud is the length of up and down, Llr is the width of left and right), wire frame center coordinates (X0, Y0), observation point coordinates (x, y), transmission current, and observation time channel.

[0033] The parameters of the layered earth model in the present embodiment are as follows: uniform half-space geoelectric structure ( p 1=100 Ω·m Typical three-layer geoelectric structure ( p 1=100 Ω·m, h 1=100 m p 2=10 Ω·m, h2 = 200 m; p 3 = 200 Ω·m Complex five-layer geoelectric structure p 1 = 100 Ω·m h 1 = 100 m p 2 = 10 Ω·m h 2 = 100 m p 3 = 100 Ω·m h 3 = 200 m p 4 = 10 Ω·m h 4 = 300 m p 5 = 200 Ω·m Fine nine-layer geoelectric structure p 1 = 100 Ω·m h 1 = 100 m p 2 = 10 Ω·m h 2 = 100 m p 3 = 100 Ω·m h 3 = 50 m p 4 = 10 Ω·m h 4 = 50 m p 5 = 200 Ω·m h 5 = 50 m p 6 = 500 Ω·m h 6 = 100 m p 7 = 300 Ω·m h 7 = 100 m p 8 = 100 Ω·m h 8 = 200 m p 9 = 200 Ω·m

[0034] The loop size is 500 m x 500 m, and the emission current is 20 A; the observation point is located at the center of the loop, that is, the line frame center coordinates and the observation point coordinates are the same; the observation time channel is logspace(log10(5*10^-5,-1,61), that is, 61 logarithmic equal interval time channels from 5*10^-5 seconds to 0.1 seconds.

[0035] If a parameter is missing in the actual situation, the default value is automatically assigned; then the persistent variable technology (persistent) is used to pre-load the Hankel transform coefficient and the sine and cosine filter coefficient, avoiding repeated loading and transmission, reducing file I / O overhead, and transmitting the above parameters and coefficients to the GPU memory, thereby reducing CPU-GPU communication transmission; at the same time, the GPU memory space is pre-allocated, the memory layout is optimized, and the gpuArray function is used to transfer the calculation data to the GPU memory.

[0036] Step two, large loop source discretization and parallel task planning: the large loop source is divided into north, south, east and west four edges, each edge is discretized into N electric dipoles (integration points) by Gauss-Legendre quadrature method; and a four-level parallel computing architecture is adopted to calculate the four edges, which fully taps the multi-core computing potential of GPU, specifically: A, loop edge parallel computing layer: the four edges (north, south, east, west) of the large loop source are allocated to different processing units (streaming multiprocessors) of GPU to realize macro task level parallel computing of each edge; the calculation of each loop edge is completely independent, and efficient load balancing is realized through parallel loop.

[0037] B, Gauss integration point parallel layer: in the calculation of each loop edge, Gauss-Legendre quadrature method is used to discretize each edge into N electric dipoles, and this calculation is performed in parallel at the GPU thread block level. This approach can fully utilize the advantages of the SIMT architecture of GPU.

[0038] C, Hankel transform parallel layer: in the calculation of the frequency domain response of a single electric dipole, the calculation of 47 Hankel transform sampling points is allocated to GPU thread block level parallel processing; through the optimized memory access mode, the efficiency of high frequency data access is ensured. The frequency domain response of a single electric dipole is calculated using the following horizontal electric dipole frequency domain magnetic field expression on the surface of a layered earth: (1) where, is the current in the electric dipole, is the length of the electric dipole, is the reflection coefficient, , is the first order Bessel function.

[0039] D, time channel parallel computing layer: for frequency-time domain conversion calculation, the calculation tasks of different time channels are parallelized to realize fine-grained data level parallel processing; through parallel processing of time channels, the overall calculation time is significantly reduced.

[0040] Step three, calculate the frequency-time domain response of each electric dipole: use the frequency-wave number-time parallel computing architecture in GPU to calculate each electric dipole in step two, so as to obtain the frequency-time domain response of each electric dipole, specifically: I, frequency domain and wave number domain sampling: log-uniform sampling is performed in the complex frequency plane and wave number domain respectively to obtain 87 frequency sampling points ωi and 47 wave number points λj.

[0041] II. Three-dimensional parameter array construction: The parameters, frequencies, and wave numbers of the layered earth model are extended to a three-dimensional parameter array (number of layers x number of wave points x number of frequency points), which is transmitted to the GPU; the purpose of this step is to reconstruct the entire calculation task into a highly regular data parallel problem.

[0042] III. Wave impedance recursive calculation reconstruction: The traditional layer-by-layer recursive calculation is converted into matrix multiplication operation to calculate the dispersion coefficient F(λ,ω) of the surface wave impedance; and the matrix calculation unit of the GPU is used to accelerate the operation, that is, each thread block or thread of the GPU independently calculates the complete recursive process under a specific (λ,ω) combination, so as to convert the serial recursion into independent parallel recursion tasks; and the bsxfun function is used instead of the repmat function for matrix operation to reduce the use of memory.

[0043] IV. Optimization of Hankel transform: The dispersion coefficient F(λ,ω) calculated in step III is used to perform Hankel transform through linear filtering algorithm and vectorized calculation method; in this process, 47-point Hankel transform coefficients are used instead of the traditional 140-point Hankel transform coefficients, so as to reduce the calculation amount under the premise of ensuring accuracy, and finally obtain the frequency domain electromagnetic field response H z (ω); V. Improved sine and cosine transform: A non-uniform sampling strategy is adopted to maintain high-density sampling in the key frequency band; the traditional 250-point sine and cosine filter coefficients are optimized to 87 points to reduce GPU-CPU data transmission; and the cubic spline interpolation algorithm is used to compensate for the accuracy loss caused by the reduction of coefficients, so as to convert the frequency domain response H z (ω) into time domain response h z (t) through sine and cosine transform; the specific formula for converting the frequency domain response H z (ω) to the time domain is as follows:

[0044]

[0045] Using digital filtering, the integral is converted to summation: (3) (4) wherein and are sine and cosine filter coefficients, is sampling interval, n is frequency sampling point index value, and k is time channel; in the above formula, that is, the time domain response h z (t).

[0046] The above process, in each task, the wave number, frequency, time channel calculation is executed according to three-layer parallel architecture.

[0047] Step four, synthesizing the large loop source response and output: after the calculation of all electric dipoles in step three is completed, the transient response of all electric dipoles is calculated and summed (Gaussian integral) on the GPU to obtain the total field response of the large loop source at the current observation point, and finally the response is transmitted from the GPU memory back to the host memory, and then the final time domain electromagnetic field value is output.

[0048] Effect verification: In order to verify the effect of the embodiment, in the forward process of the transient electromagnetic data at the current observation point in the embodiment, the existing traditional CPU serial method is used to process the transient electromagnetic data at the same observation point, and the specific results are shown in Table 1: Table 1: Comparison of calculation efficiency of the method of the application (GPU acceleration) and the traditional CPU serial method under the same conditions:

[0049] According to the comparative analysis of Table 1 and Figure 5 , 6 It can be found that the GPU accelerated forward method used in the embodiment of the application has significant advantages in calculation efficiency and accuracy compared with the traditional CPU serial algorithm: in the uniform half-space geoelectric model, the method only needs about 0.025 seconds to complete one forward simulation, which realizes a speedup of 156.9 times compared with the traditional CPU serial algorithm of 3.923 seconds; and as Figure 6 shown, the numerical solution is highly consistent with the analytical solution, and the relative errors of the early and late time channels are controlled within 4% and 0.5% respectively, which fully meets the demand of high-precision exploration. Further, Figure 7 The comparative results of the four geoelectric models in reveal a key law: as the complexity of the geologic model increases (from uniform half-space to multi-layer geoelectric structure), the parallel computing advantage of the application becomes more prominent: the speedup ratio under the nine-layer model jumps to 1213.9 times, and the GPU calculation time only increases to 0.067 seconds, while the traditional CPU time increases to 81.336 seconds, which shows that the method has excellent scalability when dealing with complex geological structures, and truly realizes the breakthrough of "the more complex the model, the more significant the efficiency improvement", providing reliable technical support for large-scale high-precision electromagnetic exploration.

[0050] The above only describes the preferred embodiments of the application, and it should be noted that for ordinary skilled persons in the art, without departing from the principles of the application, a number of improvements and refinements can be made, which should also be considered within the protection scope of the application.

Claims

1. A GPU acceleration-based large loop source transient electromagnetic fast forward method, characterized in that, The method comprises the following steps: Step one, parameter setting and pretreatment: firstly, the number of required setting parameters is determined, and each parameter is valued according to the actual situation of the current observation point, if a parameter is missing in the actual situation, a default value is automatically assigned; then the Hankel transform coefficient and the sine and cosine filter coefficient are preloaded by using the persistent variable technology, and the above parameters and coefficients are transmitted to the GPU memory; Step two, large loop source discretization and parallel task planning: the large loop source is divided into north, south, east and west four sides, each side is discretized into N electric dipoles by using the Gauss-Legendre quadrature method; and a four-level parallel computing architecture is used to calculate the four sides; Step three, calculating the frequency-time domain response of each electric dipole: in the GPU, the frequency-wavenumber-time parallel computing architecture is used to calculate each electric dipole in step two, so as to obtain the frequency-time domain response of each electric dipole; Step four, synthesizing the large loop source response and outputting: after the calculation of all electric dipoles in step three is completed, the transient response of all electric dipoles is calculated and summed on the GPU to obtain the total field response of the large loop source at the current observation point.

2. The GPU-accelerated fast forward method for transient electromagnetic method based on large loop source according to claim 1, characterized in that, The setting parameters in step one include the parameters of the layered earth model, the loop size, the wire frame center coordinates, the measurement point coordinates, the transmission current, and the observation time channel.

3. The GPU-accelerated fast forward method for transient electromagnetic method based on large loop source according to claim 1, characterized in that, The parallel task planning in step two is specifically: A, loop side parallel calculation layer: the calculation tasks of the four sides of the large loop source are distributed to different processing units of the GPU, realizing the parallel calculation of each side at the macro task level; B, Gauss integral point parallel layer: in the calculation of each loop side, the Gauss-Legendre quadrature method is used to discretize each side into N electric dipoles, and the calculation is performed in parallel at the GPU thread block level; C, Hankel transform parallel layer: in the frequency domain response calculation of a single electric dipole, the calculation of 47 Hankel transform sampling points is distributed to the GPU thread block level parallel processing; D, time channel parallel calculation layer: for frequency-time domain conversion calculation, the calculation tasks of different time channels are parallelized to realize fine-grained data level parallel processing.

4. The GPU-accelerated fast forward method for transient electromagnetic method based on large loop source according to claim 2, characterized in that, Step three is specifically: Ⅰ, frequency domain and wave number domain sampling: logarithmic uniform sampling is performed in the complex frequency plane and the wave number domain to obtain multiple frequency sampling points ωi and wave number points λj; Ⅱ, three-dimensional parameter array construction: the parameters of the layered earth model, the frequency and the wave number are expanded into a three-dimensional parameter array, and are transmitted to the GPU; Ⅲ, wave impedance recursive calculation reconstruction: the traditional layer-by-layer recursive calculation is converted into matrix multiplication operation, the three-dimensional parameter array is combined to calculate the dispersion coefficient F(λ,ω) of the ground wave impedance; and the matrix calculation unit of the GPU is used to speed up the operation, so as to convert the serial recursion into independent parallel recursion tasks; IV. Optimization of Hankel transform: using the dispersion coefficient F(λ, ω) calculated in step III, the Hankel transform is performed by a linear filtering algorithm and vectorized calculation to obtain the frequency domain electromagnetic field response H z (ω); V. Improved Sine-cosine Transform: Using non-uniform sampling strategy, high density sampling is maintained in key frequency bands, the traditional 250-point sine-cosine filter coefficient is optimized to 87 points, and the precision loss caused by coefficient reduction is compensated by combining cubic spline interpolation algorithm, so that the frequency domain response H z (ω) is converted into time domain response h z (t) by sine-cosine transform.

5. The GPU-accelerated fast forward method for transient electromagnetic method based on large loop source according to claim 4, characterized in that, In the step IV, 47-point Hankel transform coefficient is used in the Hankel transform process.

6. The GPU-accelerated fast forward method for transient electromagnetic method based on large loop source according to claim 1, characterized in that, After the step four obtains the total field response of the large loop source at the current detection point, the response is transmitted from the GPU memory back to the host memory, and then the final time domain electromagnetic field value is output.

Images