A method of reverse time migration imaging for UAV ground penetrating radar

Abstract

The invention relates to the field of electromagnetic simulation and calculation, and discloses a reverse time migration imaging method for a ground penetrating radar of an unmanned aerial vehicle, comprising: performing denoising processing on data collected by a unmanned aerial vehicle, removing interference of direct waves and reflected waves, obtaining first processed data, and constructing an initial velocity model; decoupling the Maxwell equations, and processing reflected waves by using CPML absorption boundary conditions; selecting a Yee grid to efficiently discretize a calculation area; applying a finite difference time domain method to simulate electromagnetic wave propagation, and applying a reverse time migration algorithm to perform reverse extrapolation until the electromagnetic wave field returns to an initial moment. In this process, the position of a reflection interface is calculated according to the propagation time and velocity model of the electromagnetic wave, and an underground imaging profile is obtained; finally, the imaging profiles of all measuring points are superimposed to obtain a final imaging result. The invention successfully solves the problem that a traditional imaging method based on ray migration theory in unmanned aerial vehicle GPR exploration cannot accurately image under complex surface conditions.

Description

Technical Field

[0001] The invention relates to the field of electromagnetic simulation and calculation, and in particular to a reverse time migration imaging method for an unmanned aerial vehicle ground penetrating radar. Background Art

[0002] As a non-destructive testing technology, near-ground penetrating radar has been widely used in fields such as tunnel and urban road quality testing. However, traditional ground-penetrating radar technology requires ensuring good coupling between the transmitting antenna and the ground, so its application is limited in environments such as mountainous areas with complex terrain.

[0003] For UAV ground penetrating radar signals, there are usually various surface clutters in complex surface environments, such as vegetation, buildings, roads, rocks, etc. These surface interferences will interfere with the propagation and reception of radar signals, making it difficult to detect underground targets. The strong reflection and scattering of surface clutter will confuse the target signal and reduce the resolution and accuracy of imaging. Secondly, the underground medium in a complex surface environment is usually non-uniform, including changes in soil type, water content, density, electromagnetic properties, etc. These non-uniformities will cause the radar signal to be scattered, attenuated, and refracted during propagation, causing interference to the imaging results and causing errors. Summary of the invention

[0004] The present invention aims to provide a UAV ground penetrating radar reverse time migration imaging method to effectively improve the detection and imaging accuracy of UAV GPR for shallow surface geological structures under complex surface conditions, thereby providing more accurate data support for geological structure analysis.

[0005] In order to achieve the above object, the present invention provides the following method:

[0006] The present invention provides a reverse time migration imaging method for an unmanned aerial vehicle ground penetrating radar, which is:

[0007] S1: performing a direct wave removal processing process and a data preprocessing process on the data collected by the UAV ground penetrating radar to obtain first processed data;

[0008] S2: constructing an initial velocity model based on the first processed data and geological information and prior knowledge;

[0009] S3: Decouple the Maxwell equations to adapt to the characteristics of GPR data, and use CPML absorbing boundary conditions to process the reflected waves at the boundaries to reduce the errors introduced by boundary reflections;

[0010] S4: Yee grid is selected to discretize the calculation area efficiently, and the finite difference time domain method is applied to simulate the propagation of electromagnetic waves in underground media, and the UAV GPR signals collected at each measuring point are collected to ensure stability and accuracy at various frequencies;

[0011] S5: The UAV GPR signal collected at each measuring point is reversely extrapolated using the reverse time migration algorithm until the electromagnetic wave field returns to the initial time. In this process, the position of the reflection interface is calculated according to the propagation time and velocity model of the electromagnetic wave to obtain the underground imaging profile of each measuring point. Finally, the imaging profiles of all measuring points are superimposed, integrated and accumulated to obtain the final imaging result.

[0012] Preferably, the formula for decoupling the Maxwell equations includes:

[0013]

[0014] Among them, H x and H z are the magnetic field components in the x and z directions, E y is the electric field strength in the y direction, ε is the dielectric constant, σ is the conductivity, μ is the magnetic permeability, and t is time.

[0015] Preferably, the step of introducing a CPML absorbing boundary condition to absorb the reflected wave at the boundary comprises: in a lossy medium in a two-dimensional TM mode, the general form of the electromagnetic wave frequency domain equation based on CPML is:

[0016]

[0017] Among them, si is the coordinate scaling factor.

[0018] Preferably, the conditions for selecting the Yee grid and stable high-frequency electromagnetic waves include:

[0019] Select a square Yee grid and let the high-frequency electromagnetic wave propagate along the diagonal direction of the square Yee grid. There are more than 10 sampling points in one high-frequency electromagnetic wavelength:

[0020]

[0021] Selection of the time step of the high-frequency electromagnetic wave:

[0022]

[0023] Selection of the total length of the high-frequency electromagnetic wave time:

[0024]

[0025] Where λ is the wavelength, vmax 、v min are the maximum and minimum values ​​of the electromagnetic wave propagation speed in the two-dimensional model, Δx and Δz are the vertical and horizontal grid spacings, and dt is the time step.

[0026] Preferably, the step of performing reverse time migration imaging calculation by the time-space domain iterative equation, reverse time extrapolating the GPR signal collected at each measuring point of the UAV ground penetrating radar, obtaining the imaging profile of each measuring point when the electromagnetic wave field is pushed to zero, and integrating and accumulating the imaging profile of each measuring point to obtain the final imaging comprises:

[0027] The propagation velocity v and attenuation coefficient α of electromagnetic waves can be approximately equal to:

[0028]

[0029] Under the condition that the model velocity is halved, the collected GPR signal is extrapolated in reverse time to time zero to obtain the imaging section. The model velocity should be half of the true velocity value:

[0030] v′=v / 2

[0031] It can be clearly seen from formula (1) that to ensure v′=v / 2, the relative dielectric constant needs to be replaced by ε′=4ε. At the same time, in order to make the attenuation of electromagnetic waves in unidirectional propagation consistent with the attenuation of bidirectional propagation, it can be further derived from formula (2):

[0032]

[0033] The right side of the equation is the attenuation of the two-way propagation of electromagnetic waves, while the left side is the attenuation of the one-way propagation, and ε′=4ε, which leads to:

[0034]

[0035] Substitute ε′ and σ′ into the Maxwell equations for decoupling formula to obtain the modified Maxwell equations based on reverse time migration imaging, calculate the imaging profile of each measuring point, and then integrate and accumulate the imaging profiles of each measuring point to obtain the final imaging.

[0036] The beneficial effects of the present invention are as follows: the present invention successfully applies the reverse time migration algorithm based on GPU parallel acceleration to the imaging field of unmanned aerial vehicle ground penetrating radar, and uses the algorithm to perform imaging calculations on synthetic data of shallow surface hidden danger bodies. Compared with the traditional Kirchhoff migration imaging algorithm, the algorithm can more accurately characterize the geological structure of shallow surface landslide hidden danger bodies in complex environments, and can achieve accurate imaging of shallow surface hidden danger bodies under complex surface conditions, which is of great value for early warning, prevention and control, and risk management of non-explicit landslides; through the CPML absorption boundary, compared with the split perfectly matched layer PML, CPML has a better absorption effect on grazing waves with a smaller inclination angle. When the grid size is greater than 1024, the GPU parallel algorithm has a speedup ratio of more than 10 times that of the CPU serial calculation. For large-scale FDTD calculations, the use of GPU parallel accelerated calculations can greatly improve the simulation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for the specific embodiments or the description of the prior art. In all the drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn according to the actual scale.

[0038] Figure 1 It is a wave field snapshot diagram of electric field and magnetic field components at different times provided by an embodiment of the present invention.

[0039] Figure 2 It is a flow chart of a reverse time migration imaging method for an unmanned aerial vehicle ground penetrating radar provided in an embodiment of the present invention.

[0040] Figure 3 This is a basic principle diagram of explosion reflection interface imaging of a UAV ground penetrating radar reverse time migration imaging method provided by an embodiment of the present invention.

[0041] Figure 4 This is a CPU and GPU time consumption comparison table of the embodiments of the present invention.

[0042] Figure 5 This is a comparison chart of GPU parallel and CPU serial computing efficiency according to an embodiment of the present invention.

[0043] Figure 6 It is a multi-layer geological structure model. The dark blue part is the air layer, the light purple part is the soil layer, and the reddish brown part is the bedrock layer. The yellow line segment is the actual flight trajectory, and the white line segment is the horizontal reference plane.

[0044] Figure 7 Bscan synthetic data: (a) Bscan simulation data without height correction; (b) Bscan simulation data after height correction.

[0045] Figure 8 These are the reverse time migration imaging results after velocity model and different altitude correction processing methods (a) velocity model; (b) the result of the migration calculation of Bscan simulation data without altitude correction on the horizontal reference plane; (c) the result of the migration calculation of Bscan simulation data after altitude correction on the horizontal reference plane (d) the result of the migration calculation of Bscan simulation data without altitude correction on the flight trajectory.

[0046] Fig. 9 The shallow hidden danger landslide model and the reverse time migration imaging results are shown in Figure 1. The surface is a sloping and undulating interface. (a) Relative dielectric constant model, where the yellow dots are the radar single-transmit and single-receive positions; (b) Conductivity model; (c) Synthesized Bscan data; (d) Reverse time migration imaging results. DETAILED DESCRIPTION

[0047] In order to make the technical personnel in the technical field better understand the scheme of the present invention, the technical scheme in the embodiment of the present invention will be clearly and completely described below in combination with the drawings in the embodiment of the present invention. Obviously, the described embodiment is only a part of the embodiment of the present invention, not all the embodiments. Based on the embodiment of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0048] The terms "first", "second", etc. in the specification and claims of the present invention and the above-mentioned drawings are used to distinguish different objects rather than to describe a specific order. In addition, the terms "including" and "having" and any variations thereof are intended to cover non-exclusive inclusions. For example, a process, method, device, product or end including a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units that are not listed, or may optionally include other steps or units inherent to these processes, methods, products or ends.

[0049] Reference to "embodiments" herein means that a particular feature, structure, or characteristic described in conjunction with the embodiments may be included in at least one embodiment of the present invention. The appearance of the phrase in various places in the specification does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment that is mutually exclusive with other embodiments. It is explicitly and implicitly understood by those skilled in the art that the embodiments described herein may be combined with other embodiments.

[0050] For UAV ground penetrating radar signals, there are usually various surface clutters in complex surface environments, such as vegetation, buildings, roads, rocks, etc. These surface interferences will interfere with the propagation and reception of radar signals, making it difficult to detect underground targets. The strong reflection and scattering of surface clutter will confuse the target signal and reduce the resolution and accuracy of imaging. Secondly, the underground medium in a complex surface environment is usually non-uniform, including changes in soil type, water content, density, electromagnetic properties, etc. These non-uniformities will cause the radar signal to be scattered, attenuated, and refracted during propagation, causing interference to the imaging results and causing errors.

[0051] The present invention aims to provide a reverse time migration imaging method for UAV ground penetrating radar to solve the problems of inaccurate depiction of the geological structure of shallow underground anomalies in complex environments by traditional Kirchhoff migration imaging algorithm and inaccurate imaging of shallow hidden dangers under complex surface conditions.

[0052] Embodiment 1

[0053] A specific embodiment of the present invention provides a method for reverse time migration imaging of a ground penetrating radar of an unmanned aerial vehicle, the method being as follows Figure 1 As shown, the following steps are included:

[0054] S1: performing a direct wave removal process and a data preprocessing process on the data collected by the UAV ground penetrating radar to obtain first processed data.

[0055] In the embodiment of the present invention, the mean method is used in the process of removing direct waves, and each pixel value in the image is subtracted from the mean value of the row where the pixel value is located. The main contents of data preprocessing include: data cleaning, data integration, data transformation and data reduction.

[0056] S2: constructing an initial velocity model based on the first processed data and geological information and prior knowledge.

[0057] In the embodiment of the present invention, in the reverse time migration imaging calculation of the UAV ground penetrating radar, it is necessary to pay attention to the large difference between the velocity values ​​of the air and the underground medium. Direct reverse time migration imaging calculation will cause a large number of illusions. The velocity model needs to be smoothed to reduce false reflections; under complex geological conditions, multiple reflection waves and scattered waves may cause artifacts and interference in the imaging, affecting the accuracy of the interpretation. Therefore, it is necessary to minimize multiple reflection waves and scattered waves.

[0058] S3: The Maxwell equations are decoupled to adapt to the characteristics of GPR data, and the CPML absorbing boundary conditions are used to process the reflected waves at the boundaries to reduce the errors introduced by boundary reflections.

[0059] In an embodiment of the present invention, a formula for decoupling the Maxwell equations based on a high-order finite difference method accelerated by GPU parallelization includes:

[0060]

[0061] Among them, H x and H z are the magnetic field components in the x and z directions, E y is the electric field strength in the y direction, ε is the dielectric constant, σ is the conductivity, μ is the magnetic permeability, and t is time.

[0062] The step of introducing the CPML absorption boundary condition to absorb the reflected wave at the boundary comprises:

[0063] In a lossy medium under two-dimensional TM mode, the general form of the electromagnetic wave frequency domain equation based on CPML is:

[0064]

[0065] Among them, s i is the coordinate scaling factor:

[0066]

[0067] S4: The Yee grid is selected to discretize the calculation area efficiently, and the finite difference time domain method is applied to simulate the propagation of electromagnetic waves in the underground medium. The UAV GPR signals collected at each measuring point are collected to ensure stability and accuracy at various frequencies.

[0068] In the embodiment of the present invention, the conditions for selecting the Yee grid and the stable high-frequency electromagnetic wave include:

[0069] Select a square Yee grid and let the high-frequency electromagnetic wave propagate along the diagonal direction of the square Yee grid, then the following relationship holds:

[0070] There are more than 10 sampling points in one high frequency electromagnetic wavelength

[0071]

[0072]

[0073] Selection of the time step of the high-frequency electromagnetic wave:

[0074]

[0075] Selection of the total length of the high-frequency electromagnetic wave time:

[0076]

[0077] Where λ is the wavelength, v max 、v min are the maximum and minimum values ​​of the electromagnetic wave propagation speed in the two-dimensional model, Δx and Δz are the vertical and horizontal grid spacings, and dt is the time step.

[0078] S5: The UAV GPR signal collected at each measuring point is reversely extrapolated using the reverse time migration algorithm until the electromagnetic wave field returns to the initial time. In this process, the position of the reflection interface is calculated according to the propagation time and velocity model of the electromagnetic wave to obtain the underground imaging profile of each measuring point. Finally, the imaging profiles of all measuring points are superimposed, integrated and accumulated to obtain the final imaging result.

[0079] In the embodiment of the present invention, Figure 3 As shown, the propagation velocity v and attenuation coefficient α of electromagnetic waves can be approximately equal to:

[0080]

[0081] Under the condition that the model velocity is halved, the collected GPR signal is extrapolated in reverse time to time zero to obtain the imaging section. The model velocity should be half of the true velocity value:

[0082] v′=v / 2

[0083] It can be clearly seen from formula (1) that to ensure v′=v / 2, the relative dielectric constant needs to be replaced by ε′=4ε. At the same time, in order to make the attenuation of electromagnetic waves in unidirectional propagation consistent with the attenuation of bidirectional propagation, it can be further derived from formula (2):

[0084]

[0085] The right side of the equation is the attenuation of the two-way propagation of electromagnetic waves, while the left side is the attenuation of the one-way propagation, and ε′=4ε, which leads to:

[0086]

[0087] Substitute ε′ and σ′ into the Maxwell equations for decoupling formula to obtain the modified Maxwell equations based on reverse time migration imaging, calculate the imaging profile of each measuring point, and then integrate and accumulate the imaging profiles of each measuring point to obtain the final imaging.

[0088] Embodiment 2

[0089] like Figure 4 and Figure 5As shown in the figure, the finite-difference time-domain method has a natural parallel characteristic in space: the value of the electric field (or magnetic field) at any point in space is only related to its own field value at the last moment and the field value of the magnetic field (or electric field) at the surrounding grid points in the last half time step, but has no direct relationship with other field quantities in the calculation area. Therefore, GPUs with more computing cores are used for spatial parallel processing of FDTD algorithms.

[0090] In order to test the effectiveness of the GPU-based parallel accelerated forward algorithm, the uniform air medium is divided into models with different grid sizes, and the time consumption is compared with that of the traditional CPU serial forward algorithm. Figure 4 The time consumption of different algorithms for forward modeling of models with various grid sizes at 1000 steps is shown. Figure 5 This is a comparison chart of the time consumption for forward simulation.

[0091] Depend on Figure 4 and Figure 5 It can be seen that for small grid sizes of 128 or less, the CPU is actually faster. However, when the grid size is greater than 1024, the GPU parallel algorithm has a speedup ratio of more than 10 times compared to the CPU serial calculation, indicating that the use of GPU parallel acceleration calculation for large-scale two-dimensional FDTD calculations can greatly improve the simulation efficiency.

[0092] Embodiment 3

[0093] like Figure 6 , Figure 7 As shown, Figure 6 It is a multi-layer geological structure model. The dark blue part is the air layer, the light purple part is the soil layer, and the reddish brown part is the bedrock layer. The yellow line segment is the actual flight trajectory, and the white line segment is the horizontal reference plane. Figure 7 Bscan synthetic data: (a) Bscan simulation data without height correction; (b) Bscan simulation data after height correction. Figure 6 It is a multi-layer geological structure model with a horizontal distance of 60m and a vertical distance of 25m. The model consists of three parts: the air layer, whose relative dielectric constant is 1 and the conductivity is 1mS / m; the soil layer, whose relative dielectric constant is 4 and the conductivity is 1mS / m, and the surface is a horizontal interface; the bedrock layer, whose relative dielectric constant is 8 and the conductivity is 1mS / m, and the interface between it and the soil layer is an undulating interface. The model area is divided into 600×250 square grid units, and the vertical and horizontal grid spacing is 0.1m×0.1m. In order to eliminate artificial boundary interference, the number of CPML absorption layers is set to 20 and the thickness is 2m.

[0094] The single-transmit and single-receive data acquisition mode of the ground penetrating radar equipped with the UAV is simulated, and the Ricker wavelet with a center frequency of 50MHz is used as the excitation source. The recording time is 200ns and the sampling interval is 0.167ns. Figure 6The yellow line in the middle is the actual flight trajectory of the drone obtained by GPS, which is about the horizontal distance x The function h real (x), take h0=max{h real (x)} as the horizontal reference height. Then the time correction Δτ is:

[0095]

[0096] The simulation data of this model is corrected using the flight altitude ( Figure 7 b) Processing, the result after flight altitude correction is as follows Figure 7 As shown in b. Figure 7 It can be seen from a that the phase axis of the reflected wave from the ground is curved, which is inconsistent with the horizontal surface in the model, indicating that the radar profile directly collected cannot intuitively reflect the reflected wave shape of the ground and the target below it. Figure 7 b is the radar image after flight altitude correction. By observing the corrected radar image, we can see that the direct wave is consistent with the flight trajectory, and the shape of the reflected wave from the ground is basically consistent with the undulating surface shape.

[0097] Embodiment 4

[0098] like Figure 8 As shown, Figure 8 The reverse time migration imaging results after velocity model and different height correction processing methods (a) velocity model; (b) the result of the migration calculation of Bscan simulation data without height correction on the horizontal reference plane; (c) the result of the migration calculation of Bscan simulation data after height correction on the horizontal reference plane (d) the result of the migration calculation of Bscan simulation data without height correction on the flight trajectory. This study designed three different data processing methods: (1) reverse time migration imaging calculation of Bscan simulation data without height correction on the horizontal reference plane, and then perform height correction on the imaging results; (2) migration imaging calculation of Bscan simulation data after height correction on the horizontal reference plane; (3) migration imaging calculation of Bscan simulation data without height correction on the flight trajectory. In order to reduce the false reflection waves from the ground in the reverse time migration imaging calculation, Gaussian smoothing is used to smooth the velocity model. The velocity model used is as follows: Figure 8 As shown in a, the electromagnetic wave speed of the air is set to 0.3m / ns and the speed value of the underground medium is 0.15m / ns.

[0099] The migration imaging results of the above three different methods are as follows Figure 8 The offset result of the first data processing method is shown in Figure 8As shown in b, the result of the offset imaging of the interface between the soil layer and the bedrock layer is quite different from the actual model, and the image is poorly focused in the area with large terrain fluctuations at a horizontal distance of about 40m; the second data processing method ( Figure 8 The imaging effect of c) is improved compared with the first processing method, and it can roughly outline the relatively flat area, but the reflected wave energy focusing effect is not ideal in the area with large terrain fluctuations at a horizontal distance of about 40m. Compared with the first two data processing methods, the third data processing method directly performs reverse time migration imaging from the measuring point position ( Figure 8 d) It can correctly focus the image and place the reflection boundary at the correct depth, and can still maintain a good migration effect in areas with large terrain fluctuations. Therefore, reverse time migration imaging directly from the acquisition location can obtain more accurate migration results.

[0100] Embodiment 5

[0101] like Fig. 9 As shown, Fig. 9 It is a shallow hidden danger landslide model and reverse time migration imaging results. The surface is an inclined and undulating interface. (a) Relative dielectric constant model, where the yellow line segment is the radar single-transmit and single-receive position; (b) Conductivity model; (c) Synthetic Bscan data; (d) Initial velocity model; (e) Traditional Kirchhoff migration imaging result; (e) Reverse time migration imaging result. The imaging effect of the hidden danger landslide body under the undulating terrain, a two-dimensional model of the inclined and undulating surface was constructed, with a horizontal distance of 300m and a vertical distance of 80m. The model consists of an air layer and multiple layers of water-containing media. The relative dielectric constant model ( Fig. 9 The meanings of the colors in a) are as follows: the dark blue layer is the air layer, with a relative dielectric constant of 1 and a conductivity of 0 S / m; the green layer is the gravel layer, with a relative dielectric constant of 20 and a conductivity of 0.06 S / m; the light blue layer is the clay layer, with a relative dielectric constant of 12 and a conductivity of 0.05 S / m; the dark yellow layer is the sandstone layer, with a relative dielectric constant of 25 and a conductivity of 0.01 S / m; the light yellow layer is the mudstone layer, which has a large internal water content, a relative dielectric constant of 35 and a conductivity of 0.01 S / m.

[0102] The relative permittivity model and conductivity model of the model are as follows: Fig. 9 (a) and Fig. 9As shown in (b), in order to simulate the real landslide hazard, the surface is set to an inclined and undulating terrain. The multiple geological bodies in the model are all water-containing soils. When the water content increases, the area will become loose, which will in turn cause landslides. The grid division, grid spacing, and CPML layer number settings are consistent with the previous model. The simulated drone is equipped with a ground penetrating radar with a single transmission and single reception data acquisition method. The Ricker wavelet with a main frequency of 50MHz is also used as the excitation source, with a recording time of 600ns and a sampling interval of 0.167ns. The flight acquisition trajectory is shown in Fig. 9 (a) Shown by the yellow line segment.

[0103] Simulation Bscan synthetic data such as Fig. 9 As shown in c, due to the absorption of electromagnetic waves by the water body, the electromagnetic reflection wave from the lower interface of the water body is severely attenuated. The proposed reverse time migration imaging calculation is directly performed on the flight trajectory to obtain the migration result ( Fig. 9 f) and compared with the Kirchhoff migration imaging results. It can be found that the two migration results have the same surface. Fig. 9 The topographic undulation in a is basically consistent with that in the figure, and the interface between the gravel soil layer and the clay layer is Fig. 9 The two figures are basically consistent, and both figures can roughly outline the contour of the upper interface of the mudstone layer, but the red dotted circle and the arrow indicate that the reverse time migration imaging algorithm of the UAV ground penetrating radar can accurately image the underground structure of the shallow surface potential landslide body under the undulating terrain.

[0104] The above is only an embodiment of the present invention, and the common knowledge such as the known specific technical solutions or characteristics in the solution is not described in detail here. It should be pointed out that for those skilled in the art, without departing from the solution of the present invention, several modifications and improvements can be made, which should also be regarded as the protection scope of the present invention, and these will not affect the effect of the present invention and the practicality of the patent. The scope of protection required by this application shall be based on the content of its claims, and the specific implementation methods and other records in the specification can be used to interpret the content of the claims.

Claims

1. A method for reverse time migration imaging of unmanned aerial vehicle ground penetrating radar, characterized in that: The method comprises: S1: performing a direct wave removal processing process and a data preprocessing process on the data collected by the UAV ground penetrating radar to obtain first processed data; S2: constructing an initial velocity model based on the first processed data and geological information and prior knowledge; S3: Decouple the Maxwell equations to adapt to the characteristics of GPR data, and use CPML absorbing boundary conditions to process the reflected waves at the boundaries to reduce the errors introduced by boundary reflections; S4: Yee grid is selected to discretize the calculation area efficiently, and the finite difference time domain method is applied to simulate the propagation of electromagnetic waves in underground media, and the UAV GPR signals collected at each measuring point are collected to ensure stability and accuracy at various frequencies; S5: The UAV GPR signal collected at each measuring point is reversely extrapolated by a reverse time migration algorithm until the electromagnetic wave field returns to the initial time. In this process, the position of the reflection interface is calculated according to the propagation time and velocity model of the electromagnetic wave to obtain the underground imaging profile of each measuring point. Finally, the imaging profiles of all measuring points are superimposed, integrated and accumulated to obtain the final imaging result. The conditions for selecting the Yee grid include: Select a high-frequency electromagnetic wavelength with more than 10 sampling points, Selection of the time step of the high frequency electromagnetic wavelength: Selection of the total time length of the high-frequency electromagnetic wavelength: Where λ is the wavelength, Vmax and Vmin are the maximum and minimum values ​​of the electromagnetic wave propagation speed in the two-dimensional model, Δx and Δz are the vertical and horizontal grid spacings, and dt is the time step; The UAV GPR signal collected at each measuring point is reversely extrapolated by applying a reverse time migration algorithm until the electromagnetic wave field returns to the initial time. In this process, the position of the reflection interface is calculated according to the propagation time and velocity model of the electromagnetic wave to obtain the underground imaging profile of each measuring point. Finally, the imaging profiles of all measuring points are superimposed, integrated and accumulated to obtain the final imaging result, including the following steps: The propagation velocity v and attenuation coefficient α of electromagnetic waves are approximately equal to: Where ε is the dielectric constant; Under the condition that the model velocity is halved, the collected GPR signal is extrapolated in reverse time to time zero to obtain the imaging section. The model velocity should be half of the true velocity value: v′=v / 2 (10) From formula (8), we can get that to ensure v′=v / 2, the relative dielectric constant needs to be replaced by ε′=4ε. At the same time, in order to make the attenuation of electromagnetic waves in unidirectional propagation consistent with the attenuation of bidirectional propagation, we can deduce from formula (9): The right side of the equation is the attenuation of the two-way propagation of electromagnetic waves, while the left side is the attenuation of the one-way propagation, and ε′=4ε, which leads to: Substitute ε′ and σ′ into the Maxwell equations for decoupling, and obtain a modified Maxwell equation group based on reverse time migration imaging. Calculate the imaging profile of each measuring point, and then integrate and accumulate the imaging profiles of each measuring point to obtain the final imaging.

2. The method for reverse time migration of ground penetrating radar using an unmanned aerial vehicle according to claim 1, characterized in that: The formula for decoupling the Maxwell equations includes: Among them, H x and H z are the magnetic field components in the x and z directions, E y is the electric field strength in the y direction, ε is the dielectric constant, σ is the conductivity, μ is the magnetic permeability, and t is time.

3. The method for reverse time migration imaging of an unmanned aerial vehicle ground penetrating radar according to claim 1, characterized in that: The step of using CPML to absorb the reflected wave at the boundary condition includes: In a lossy medium under two-dimensional TM mode, the general form of the electromagnetic wave frequency domain equation based on CPML is: The si in formula (2) is the coordinate scaling factor: Among them, H x and H z are the magnetic field components in the x and z directions, E y is the electric field strength in the y direction, and μ is the magnetic permeability.

Images