CUDA (Compute Unified Device Architecture) and PyTorch mixed compiled ground penetrating radar rapid full waveform inversion method

Abstract

The invention is suitable for the technical field of ground penetrating radar exploration, and provides a CUDA (Compute Unified Device Architecture) and PyTorch mixed compiled ground penetrating radar rapid full waveform inversion method, a PyTorch automatic differential function is self-defined by innovatively carrying out joint compilation on PyTorch and CUDA kernel functions, and a PyTorch platform is combined with CUDA mixed programming to realize automatic updating of two parameters of dielectric constant and conductivity. According to the method, GPU parallel calculation is used, the calculation efficiency of the ground penetrating radar data full-waveform inversion technology is remarkably improved, a two-parameter model can be quickly inverted, and the overall calculation efficiency is improved. The improvement of the efficiency is of great significance to the improvement of the overall effect of ground penetrating radar exploration and the practical application in the fields of urban underground space detection, road maintenance and the like.

Description

Technical Field

[0001] The present invention belongs to the technical field of ground penetrating radar exploration, and particularly relates to a fast full-waveform inversion method for ground penetrating radar with CUDA and PyTorch hybrid compilation. Background Technique

[0002] Ground penetrating radar is an effective means for detecting underground targets developed in recent decades. It is a non-destructive detection technology. Compared with other conventional underground detection methods, it has the advantages of fast detection speed, continuous detection process, high resolution, convenient and flexible operation, low detection cost, etc., and is increasingly widely used in the field of engineering investigation.

[0003] Traditional ground penetrating radar inversion relies on travel-time inversion or ray theory. The resolution is limited by the wavelength, and centimeter-scale structures (such as pipeline cracks and rock and soil fissures) cannot be accurately identified. Moreover, the inversion error is significant in multi-media complex strata. Full-waveform inversion can provide fine underground structure information, but its low computational efficiency and long time consumption limit its practical application. First, the full-waveform inversion of ground penetrating radar needs to solve a large-scale non-linear optimization problem. The parameter space increases exponentially with the refinement of the model grid, resulting in a huge amount of calculation. Each iteration requires solving the wave equation multiple times, and the computational cost is high. Second, the complexity of numerical simulation is also one of the restrictive factors. The reason for the high cost of forward simulation is that it relies on methods such as finite difference or finite element to solve Maxwell's equations. It is restricted by the Courant-Friedrichs-Lewy condition, which is a necessary condition for the convergence of explicit time marching schemes when numerically solving partial differential equations. Its core idea is that the time step in numerical calculation must be less than the time required for physical disturbances to propagate through the grid cells, otherwise the simulation results will be distorted. Therefore, high requirements are imposed on grid discretization. Therefore, the full-waveform inversion of ground penetrating radar requires high-performance computing support, and the traditional CPU-based full-waveform inversion method for ground penetrating radar data has insufficient parallel efficiency. Summary of the Invention

[0004] The purpose of the embodiments of the present invention is to provide a fast full-waveform inversion method for ground penetrating radar with CUDA and PyTorch hybrid compilation, aiming to solve the problems proposed in the above background technique.

[0005] The embodiments of the present invention are implemented as follows. The fast full-waveform inversion method for ground penetrating radar with CUDA and PyTorch hybrid compilation includes the following steps: Step 1: Forward reconstruction; The forward propagation is from time 0 to time T. Step 1.1: First, store the receiver signals at the current time; Step 1.2: Electromagnetic wave field propagation; In the full-wave inversion of ground penetrating radar, the forward simulation uses the finite-difference time-domain method to describe the propagation behavior of electromagnetic waves in the medium by solving Maxwell's equations. The propagation of electromagnetic waves in ground penetrating radar is described by the time-domain Maxwell curl equation: Faraday's law of electromagnetic induction: ; where, is the curl operation; E is the electric field strength (V / m); is the magnetic permeability (H / m); H is the magnetic field strength (A / m); is the equivalent magnetic loss (Ω / m); t is the time (s); Ampere-Maxwell law: ; where, is; is the permittivity (F / m).

[0006] The finite-difference time-domain method discretizes Maxwell's equations in space and time using a Yee grid to obtain the update equations: ; ; ; where, is the grid step size in the X direction, is the grid step size in the Y direction, is the time step size, which needs to satisfy the Courant stability condition, and c is the speed of light in the medium; is the value of the electric field in the Z direction at time n; is the value of the magnetic field in the Y direction at time n; is the value of the magnetic field in the X direction at time n; i , j are the Yee grid coordinates; Update the magnetic field and its PML first according to the above equations, and then update the electric field and its PML; Step 1.3: Excitation source; The excitation source uses a Ricker wavelet with a main frequency of 50 MHz - 1 GHz: ; where, is the amplitude value of the waveform at time t, dimensionless; is the center frequency of the wavelet, unit: Hz; e and π are constants; Step 1.4: Store the wave fields required for inversion. Step 2: Inversion; The adjoint wavefield propagates from time T to time 0; Step 2.1: Excite the adjoint field source; Forward propagation obtains the forward modeling result and the forward modeling wavefield, and the corresponding residual is obtained using the specified loss function as the source of the adjoint wavefield for the backpropagation part. This step is implemented by PyTorch; Step 2.2: Propagate the adjoint electromagnetic wavefield; First, propagate the electric field and its PML, and then propagate the magnetic field and its PML.

[0007] Step 2.3: Calculate the gradients of the relative permittivity and conductivity through the following formulas; ; ; where, is the gradient of the relative permittivity, is the gradient of the conductivity, E z is the electric field in the z direction, is the adjoint field, N s is the number of adjoint sources, and T is the time window; Step 3: Update the model; After the forward modeling part and the inversion part are completed, the gradients are obtained and the initial model is updated using the PyTorch optimizer until the stopping condition is reached.

[0008] The fast full-waveform inversion method for ground penetrating radar with CUDA and PyTorch hybrid compilation provided by the embodiments of the present invention innovatively combines PyTorch with CUDA kernel functions and customizes the PyTorch automatic differentiation function, achieving the following remarkable effects in the full-waveform inversion of ground penetrating radar: A breakthrough improvement in computing performance. Compared with traditional CPU-based methods, the full-waveform inversion efficiency is increased by dozens of times, and the single iteration time is shortened from dozens of seconds to several seconds; Improved usability. The Python and C++ hybrid architecture balances development efficiency and computing performance, integrates into the PyTorch ecosystem, and directly calls PyTorch optimizers, loss functions, etc.

[0009] In the simulation, this method shows excellent performance. It can quickly invert the two-parameter model, improving the overall computing efficiency. This improvement in efficiency is of great significance for enhancing the overall effect of ground penetrating radar exploration and its practical applications in fields such as urban underground space detection and road maintenance. Brief Description of the Drawings

[0010] Figure 1Flowchart of the fast full-waveform inversion method for ground penetrating radar with CUDA and PyTorch hybrid compilation provided by the embodiments of the present invention; Figure 2 It is a single-parameter inversion result diagram, where, Figure 2 a is the dielectric constant ε Initial model, Figure 2 b is the dielectric constant of this method ε Prediction result; Figure 2 c is the dielectric constant ε True model; Figure 2 d is the dielectric constant prediction result of the traditional method ε Prediction result; Figure 3 It is a two-parameter inversion result diagram; where, Figure 3 a is the initial model of the relative dielectric constant ; Figure 3 b is the prediction result of the relative dielectric constant ; Figure 3 c is the true model of the relative dielectric constant ; Figure 3 d is the initial model of the conductivity σ ; Figure 3 e is the prediction result of the conductivity σ ; Figure 3 f is the true model of the conductivity σ . Specific implementation

[0011] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0012] The following describes the specific implementation of the present invention in detail with reference to specific embodiments.

[0013] torch.autograd.Function is one of the core components of the PyTorch automatic differentiation system and is used to define the custom forward propagation and backward propagation calculation logics.

[0014] As Figure 1 shown, it is a fast full-waveform inversion method for ground penetrating radar with CUDA and PyTorch hybrid compilation provided by an embodiment of the present invention. In this method, the forward modeling is reconstructed as the forward propagation part, and the backward propagation part is to solve the gradient by the adjoint state method. Specifically: Step 1: Forward modeling reconstruction; The forward modeling propagates from the 0th moment to the Tth moment; Step 1.1: First, store the receiver signal at the current moment; Step 1.2: Electromagnetic wave field propagation; In the full-wave inversion of ground penetrating radar, the forward simulation uses the finite-difference time-domain method to describe the propagation behavior of electromagnetic waves in the medium by solving Maxwell's equations. The propagation of electromagnetic waves in ground penetrating radar can be described by the time-domain Maxwell curl equation (ignoring the high-order terms of displacement current): Faraday's law of electromagnetic induction (magnetic field update): ; where is the curl operation; E is the electric field strength (V / m); is the magnetic permeability (H / m); H is the magnetic field strength (A / m); is the equivalent magnetic loss (Ω / m); t is the time (s); Ampere-Maxwell's law (electric field update): ; where is; is the permittivity (F / m).

[0015] The finite-difference time-domain method uses the Yee grid to discretize Maxwell's equations in space and time to obtain the update equations: ; ; ; where is the grid step size in the X direction, is the grid step size in the Y direction, is the time step size, which needs to satisfy the Courant stability condition, and c is the speed of light in the medium; is the value of the electric field in the Z direction at the nth moment; is the value of the magnetic field in the Y direction at the nth moment; is the value of the magnetic field in the X direction at the nth moment; i and j are the Yee grid coordinates.

[0016] At the same time, the ground penetrating radar simulation needs to process the perfectly matched layer (PML). An artificial loss layer is introduced at the boundary to absorb waves to prevent wave reflection back to the target wave field and affect the results.

[0017] Update the magnetic field and its PML first, and then update the electric field and its PML according to the above equations.

[0018] Step 1.3: Excitation source; The commonly used excitation source is the Ricker wavelet (dominant frequency 50 MHz–1 GHz): ; where, is the amplitude value of the waveform at time t, dimensionless; is the center frequency of the wavelet, unit: Hz; e and π are constants.

[0019] Step 1.4: Store the wave field required for inversion.

[0020] Step 2: Inversion; The adjoint wave field propagates from time T to time 0.

[0021] Step 2.1: Excite the adjoint field source; The forward propagation obtains the forward modeling result and the forward modeling wave field, and the corresponding residual is obtained using the specified loss function as the source of the adjoint wave field for the backward propagation part. This step is implemented by PyTorch.

[0022] Step 2.2: Propagate the adjoint electromagnetic wave field; First, propagate the electric field and its PML, and then propagate the magnetic field and its PML.

[0023] Step 2.3: Calculate the gradients of the relative permittivity and conductivity through the following formulas.

[0024] ; ; where, is the gradient of the relative permittivity, is the gradient of the conductivity, E z is the electric field in the z direction, is the adjoint field, N s is the number of adjoint sources, and T is the time window.

[0025] Step 3: Update the model; After the forward modeling part and the inversion part are completed, the gradients are obtained and the initial model is updated using the PyTorch optimizer until the stopping condition is reached.

[0026] In the embodiments of the present invention, all forward wavefield propagation, implementation of the perfectly matched layer, adjoint wavefield propagation, and gradient calculation are accelerated using CUDA kernel functions written in C++, significantly improving the calculation speed and embedding it into automatic differentiation. CUDA kernel functions are the core units of GPU parallel computing. Grid, Block, and Thread are the three-layer hierarchical structure of parallel computing, used to organize large-scale parallel computing tasks on the GPU. Their relationship is similar to a three-dimensional "task decomposition map", and its workflow can be divided into the following stages: 1) The host (CPU) calls the CPU to start the kernel function through the <<<grid, block>>> syntax, specifying the parallel hierarchy (grid + thread block); 2) Thread allocation, the GPU creates a thread grid according to gridDim and blockDim, and each thread independently executes the kernel function code; 3) Parallel computing, all threads synchronously execute the kernel function, and tasks are distinguished by the thread ID (threadIdx / blockIdx); 4) Result write-back, the calculation result is written into the global memory (Global Memory) for subsequent CPU or GPU operations.

[0027] CUDA kernel functions can achieve dozens of times of acceleration compared with the CPU solution in the full-waveform inversion of ground-penetrating radar. Their high performance stems from the deep optimization of the GPU hardware architecture and parallel computing design. The reason is that the GPU has a large number of parallel computing units. The CPU usually has 4 to 64 cores, which are suitable for serial and complex logical tasks. While the GPU has thousands of CUDA cores, which are designed specifically for parallel computing and use the single-instruction multiple-thread architecture. All CUDA cores synchronously execute the same instruction but process different data, such as calculating the electric field update of all points in the finite-difference time-domain grid simultaneously, which perfectly adapts to the electromagnetic field update of the finite-difference time-domain network.

[0028] Compared with the existing methods based on the single PyTorch platform's GPU automatic differentiation framework, the custom automatic differentiation strategy based on CUDA kernel functions adopted in the present invention shows significant advantages in multiple aspects: First, in terms of computational efficiency and performance optimization, by explicitly defining the forward and backward propagation processes, this method avoids the memory inflation problem brought by PyTorch's default automatic differentiation mechanism in solving wave equations, thus achieving efficient adaptation and in-depth optimization for the full-waveform inversion task of ground-penetrating radar; Second, the key numerical calculation modules are implemented by highly parallel CUDA kernel functions, which can fully exploit the acceleration potential of the GPU architecture in processing large-scale grid calculations and significantly improve the overall computational performance; Third, in terms of development efficiency and ecological compatibility, this method is highly integrated with the PyTorch ecosystem, supporting direct invocation of native tools such as torch.optim optimizers without the need to re-implement relevant modules, reducing development costs; Finally, relying on the hybrid programming paradigm of Python and C++, this framework realizes rapid implementation of algorithm logic at the high-level language level and extreme optimization of performance-critical paths through CUDA kernel functions at the bottom layer, combining development flexibility and runtime efficiency, providing an extensible and maintainable engineering solution for complex electromagnetic inversion problems.

[0029] As a preferred embodiment of the present invention, applying this method to the model described below achieves two-parameter inversion, while the same type of method using CPU only achieves single-parameter inversion.

[0030] Under the condition of the same grid size, it runs on the same computer equipped with an intel i7-12700k processor, 64GB of memory, and an NVIDIA RTX4090 graphics card (24GB of video memory). Compared with the same type of method using CPU, a single-parameter test was carried out on the OverThrust model with a grid size of 100*200, a spatial sampling interval of 0.02m, a time sampling interval of 0.04ns, and the source and receiver moving 0.04m to the right with a step size of 0.04m above, for a total of 100 steps.

[0031] The time spent for a single iteration is as shown in Table 1 below. The running time is 237.53s, and a total of 100 iterations are run.

[0032] Table 1 Comparison between this method and other methods ; The single-parameter inversion results are as Figure 2 shown, where Figure 2 a is the dielectric constant ε initial model, Figure 2 b is the dielectric constant of this method ε prediction result; Figure 2 c is the dielectric constant εThe true model; Figure 2 d is the dielectric constant of the traditional method ε The prediction result.

[0033] Meanwhile, two-parameter inversion was also carried out. The grid size was 200*100, the spatial sampling interval was 0.05 m, the time sampling interval was 0.1 ns, there were a total of 40 sources equally spaced on the left side of the model, and there were a total of 100 receivers equally spaced on the right side of the model. The running time was 480 s, and a total of 240 iterations were run. The traditional full waveform inversion takes dozens of minutes.

[0034] The results of the two-parameter inversion are as Figure 3 shown, where Figure 3 a is the initial model of the relative dielectric constant ; Figure 3 b is the prediction result of the relative dielectric constant ; Figure 3 c is the true model of the relative dielectric constant ; Figure 3 d is the initial model of the conductivity σ ; Figure 3 e is the prediction result of the conductivity σ ; Figure 3 f is the true model of the conductivity σ ;

[0035] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A fast full waveform inversion method for ground penetrating radar compiled by CUDA and PyTorch, characterized in that: The following steps are involved: Step 1: Forward reconstruction, simulating the propagation of electromagnetic waves in the medium through the finite-difference time-domain method; Step 2: Inversion; The accompanying wave field propagates from time T to time 0; Step 2.1: Excite the adjoint field source; forward propagate to obtain the forward modeling result and the forward modeling wave field, and use the specified loss function to obtain the corresponding residual as the source of the adjoint wave field for the backward propagation part. This step is implemented by PyTorch; Step 2.2: Accompanying the propagation of electromagnetic wave fields; first propagate the electric field and its PML, then propagate the magnetic field and its PML; Step 2.3: Calculate the dielectric constant and conductivity gradient using the following formulas; ； ； in, is the gradient of the relative dielectric constant, is the conductivity gradient, E z is the electric field in the z direction, To accompany the field, N s is the number of companion sources, T It is the time window; Step 3: Update the model; After the forward and inversion parts are completed, the gradients are obtained and the PyTorch optimizer is used to update the initial model until the stopping condition is reached.

2. The CUDA and PyTorch hybrid compiled ground penetrating radar fast full waveform inversion method according to claim 1 is characterized in that: The step 1 comprises the following specific steps: Forward modeling propagates from time 0 to time T; Step 1.1: First, store the receiver signal at the current moment; Step 1.2: Electromagnetic wave field propagation; In the full waveform inversion of ground penetrating radar, the forward simulation uses the time-domain finite difference method to describe the propagation behavior of electromagnetic waves in the medium by solving the Maxwell equations; the electromagnetic wave propagation of ground penetrating radar is described by the time-domain Maxwell curl equation: Faraday's law of electromagnetic induction: ； in, is the curl operation; E is the electric field strength, unit: V / m; is the magnetic permeability, unit: H / m; H is the magnetic field intensity, unit: A / m; is the equivalent magnetic loss, unit: Ω / m; t is time, unit: s; Ampere-Maxwell Laws: ； in, is the conductivity; is the dielectric constant, unit: F / m; The finite-difference time-domain method uses Yee grids to discretize Maxwell's equations in time and space, and obtains the updated equations: ； ； ； in, is the grid step size in the X direction, is the grid step size in the Y direction, is the time step, which must satisfy the Courant stability condition, and c is the speed of light in the medium; is the value of the electric field in the Z direction at time n; is the value of the magnetic field in the Y direction at time n; is the value of the magnetic field in the X direction at time n; i , j is the Yee grid coordinate; According to the above equation, first update the magnetic field and its PML, then update the electric field and its PML; Step 1.3: Excitation source: The excitation source uses Ricker wavelet, with a main frequency of 50 MHz-1 GHz: ； in, is the amplitude value of the waveform at time t, dimensionless; is the center frequency of the wavelet, unit: Hz; e and π is a constant; Step 1.4: Store the wave fields to be used in the inversion.

Images