A layered cross-correlation backprojection imaging method for ground penetrating radar based on ray theory

Abstract

The present invention discloses a layered cross-correlation back-projection imaging method for ground-penetrating radar based on ray theory. The method comprises the following steps: first, according to Snell's law, a layered medium ray propagation model under a common offset system is established; then, the two-way delay on the propagation model is calculated for coherent superposition in the subsequent imaging process; based on the ray propagation model and the two-way delay at the refraction point, the two-way delay information of any point on the transmission path is quickly calculated; then, in the entire area to be imaged, only the two-way delay information on the transmission path is known; in order to avoid low energy due to too few coherent superposition times of deep targets during the imaging process, the two-way delay matrix is ​​filled with a default matrix; finally, the filled delay matrix provides echo energy information during the imaging point superposition process, and the entire imaging area is subjected to back-projection coherent superposition to obtain a final high-precision imaging result.

Description

Technical Field

[0001] The present invention relates to a layered cross-correlation back-projection imaging method of a ground-penetrating radar based on ray theory. The method is applicable to ground-penetrating radar imaging processing of shallow layers and belongs to the technical field of radar signal processing. Background Art

[0002] Ground Penetrating Radar (GPR) is a specialized radar commonly used for nondestructive detection of targets within dielectric media. It transmits an ultra-wideband signal underground and, based on reflection and scattering of electromagnetic waves at discontinuities in the subsurface, can infer information such as the location and morphology of targets within the medium. GPR operates at frequencies generally between 1 MHz and 10 GHz, providing high resolution and imaging efficiency. In recent years, GPR technology has been widely used in road inspection, archaeological research, glacial exploration, military reconnaissance, urban construction, and medical imaging.

[0003] Interpreting raw GPR data is challenging, making it difficult to directly observe target information. This is why GPR imaging technology has emerged. GPR imaging technology can offset the echo energy in the raw data to its true physical location, enabling intuitive visualization of underground target information. This significantly reduces the difficulty and uncertainty of data interpretation and improves the accuracy of target detection and identification. Classic GPR imaging algorithms include back projection (BP), Kirchhoff integral migration, range migration (RM), and reverse time migration (RTM).

[0004] In most cases, GPR applications involve layered media structures rather than single media structures. For example, asphalt pavements have a surface layer, base layer, and cushion layer; building exterior walls have a plaster layer, insulation layer, and wall structure; and soil also exhibits stratification upwards due to the influence of moisture content, density, and material. Classic GPR imaging techniques are often only suitable for imaging single-media scenarios. Therefore, research on ground-penetrating radar imaging methods for layered media is of great practical significance for improving GPR imaging quality, preventing defocusing of the target area, and enhancing ground-penetrating radar detection performance. The impact of layered media on electromagnetic wave propagation is primarily manifested in variations in wave velocity due to varying dielectric constants within the layers, and discontinuities in electromagnetic waves due to the interfaces between layers.

[0005] Currently, GPR layered imaging methods are mainly divided into three categories: time-domain algorithms, frequency-domain algorithms, and electromagnetic backscattering inversion. Time-domain algorithms are mainly represented by the Layer Back Projection (LBP) imaging algorithm. This algorithm analyzes the variation of the refraction point position to solve the two-way delay at the imaging location. By introducing an approximate refraction point equation instead of the complex solution of a quartic equation, it significantly reduces the computational effort. However, this algorithm is only applicable to two-layer media and fails when the refractive index is less than 1. Frequency-domain algorithms are mainly represented by the Layer Range Migration (LRM) imaging algorithm. This algorithm derives the longitudinal wave equation and resolves the electromagnetic wave discontinuity problem at the layered media by extending the electromagnetic wave phase. It then performs range migration imaging of each layer interface. This algorithm uses fast Fourier transform technology and has a faster computational speed, but its imaging accuracy is inferior to that of time-domain algorithms. Furthermore, because the imaging results are stitched together in the range direction, the imaging quality of targets across layers is low. Electromagnetic inverse scattering algorithms can solve the scattered field of three-dimensional layered media by studying the Green's function of the layered medium, thereby inverting the target's true physical size and properties. However, they are prone to falling into local optimal solutions and are computationally intensive, making real-time imaging difficult. Therefore, research on fast, high-quality imaging of layered media by ground-penetrating radar is essential. Summary of the Invention

[0006] The purpose of the present invention is to overcome the problem that existing ground penetrating radar imaging methods are difficult to quickly and accurately image layered medium scenes. A ground penetrating radar layered cross-correlation back-projection (RLCBP) imaging method based on ray theory is proposed. The method comprises the following steps: firstly, a layered medium ray propagation model under a common offset system is established according to Snell's law; secondly, the two-way delay on the propagation model is calculated for coherent superposition in the subsequent imaging process; since the two-way delay of an imaging point in a uniform medium is positively correlated with the distance, the two-way delay information of any point on the transmission path can be quickly calculated based on the ray propagation model and the two-way delay at the refraction point; secondly, in the entire area to be imaged, only the two-way delay information on the transmission path is known, and the two-way delay information has the characteristics of being dense on the top and sparse on the bottom, dense in the middle and sparse on both sides; in order to avoid too few coherent superposition times of deep targets during the imaging process, resulting in low energy, the two-way delay matrix is ​​filled with a default matrix; finally, the filled delay matrix provides echo energy information in the imaging point superposition process, and the entire imaging area is subjected to back-projection coherent superposition, and a signal cross-correlation operation is introduced to suppress sidelobe and clutter energy, thereby obtaining a final high-precision imaging result.

[0007] The beneficial effects of the present invention are:

[0008] The present invention is applied in the field of ground penetrating radar imaging, and solves the problem of fast and accurate imaging of ground penetrating radar in layered scenes. The refraction points at the layered interface are often difficult to solve and require a large amount of calculation. Classical BP imaging reduces the amount of calculation by introducing an approximate formula for the refraction point, but it is only applicable to double-layer media. The present invention quickly calculates the two-way time delay information of any point in the layered scene by establishing a ray propagation model, avoiding the complex refraction point calculation and solution process, and at the same time realizing the imaging of multi-layer media. From the perspective of calculation time, the calculation time of time domain algorithms is often greater than that of frequency domain algorithms. The present invention greatly improves the calculation efficiency and avoids the problem of long calculation time of time domain algorithms; at the same time, it ensures that the time domain algorithm has the advantage of high-resolution imaging accuracy, introduces signal cross-correlation operations, reduces the impact of side lobes and clutter on imaging quality, and realizes accurate and fast imaging of layered media. BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 is a signal processing flow chart of an embodiment of the present invention;

[0010] Figure 2 It is a simulation experiment scene diagram in the method of the present invention;

[0011] Figure 3 This is a result diagram of the present invention imaging the simulation scene.

[0012] Specific implementation process

[0013] The purpose of this invention is to overcome the shortcomings of existing ground penetrating radar imaging algorithms, to a certain extent solve the problem of difficulty in imaging in layered media by time domain ground penetrating radar imaging algorithms, and greatly improve the calculation speed and ensure that the imaging results have the advantage of high resolution.

[0014] The present invention is achieved by the following steps:

[0015] Step 1: Establish a layered medium ray propagation model.

[0016] Ray theory is a high-frequency approximation of the wave equation. When the frequency of an electromagnetic wave is high enough, its propagation can be viewed as a ray. The ray propagation equation is also a simplified version of the wave equation. Therefore, when analyzing the propagation of electromagnetic waves in a medium from a time domain perspective, the electromagnetic wave can be viewed as propagating in the form of rays. Snell's law is used to solve the propagation path in layered media, and a ray propagation model for layered media can be established.

[0017] The receiving antenna and transmitting antenna of the ground penetrating radar are co-located and move at a fixed distance along the x-direction to measure directly downward to obtain B-scan data in the spatial and temporal dimensions. The position of the transmitting and receiving antennas is regarded as the emission source of the ray, and each transmitted signal has a certain beam width. , with a fixed angle step size The beam Decompose into ray transmission paths. The propagation model of the ray is regarded as the signal propagation path of a ground penetrating radar measurement. The antenna position is The signal propagation path It can be expressed as:

[0018]

[0019] in Representative starting point The antenna transmits Ray trajectory, Is the angle between the ray direction and the antenna measurement direction. In a layered medium scenario, the thickness and dielectric constant of each layer will affect the transmission process of the ray. The layer thickness of the layered structure is , the dielectric constant is When the intersection of the ray at the interface of the medium and the two-way delay of the signal receiving and transmitting process at this position are calculated according to Snell's law. For example, The composition of the segment ray:

[0020] ;

[0021] in and They are The starting and ending positions of the segment ray, and It is The length and angle of the segment ray, The end point of the previous sub-ray is the starting point of the next sub-ray, and the following physical relationship exists:

[0022] ;

[0023] ;

[0024] According to the layer thickness and dielectric constant, the Duanzi ray physics information Solve the Duanzi ray physics information , the following physical relationship exists:

[0025] ;

[0026] ;

[0027] Antenna position and ray initial angle It is known in the initial assumption and can be solved recursively layer by layer based on the information of the known layered scenario, thus obtaining the ray propagation model in the entire layered medium scenario.

[0028] Step 2: Generate a round-trip delay matrix based on the ray propagation model.

[0029] The basic concept of the BP algorithm is to divide the imaging area into discrete grid points and perform coherent superposition of the backscattered echoes at each imaging point to obtain the scattering intensity information at that imaging point. Coherent superposition is only possible if the round-trip delay of the imaging point relative to each synthetic aperture is known. Therefore, it is essential to obtain the round-trip delay matrix between the antenna and the imaging area. This step typically requires a considerable amount of computational time. The ray propagation model provides the positional information of signal propagation. The round-trip delay information on the rays must be calculated and applied to the imaging area to generate the round-trip delay matrix.

[0030] Knowing the layered medium scene information and using it to find the corresponding ray propagation model, the two-way delay matrix of the area to be imaged can be solved. First, solve Round-trip delay at the end point of the segment ray :

[0031]

[0032]

[0033] in Indicates that electromagnetic waves The propagation speed in the layer medium is such that the electromagnetic wave needs to go back and forth twice on the path to reach each point. Therefore, when calculating the delay between sending and receiving, the wave speed needs to be set to half of the original one. Since each layer is a uniform medium, the distance of the sub-ray in this layer is directly proportional to the round-trip delay. The time delay vector of the sub-ray path is obtained by discretizing the step size , Representing coordinates The round-trip delay at In order to avoid redundant calculations caused by too many discrete points, the step size Make this request:

[0034]

[0035] in is the wave velocity in the medium, By discretizing the sub-rays, we can find the physical position and the corresponding round-trip delay of any discrete point on the entire ray propagation path, and convert the ray propagation model into the round-trip delay matrix of the imaging grid:

[0036]

[0037] in It is a two-dimensional matrix data, indicating that the coordinates in the imaging grid matrix are Round trip delay information at , Indicates the first The first track The two-way delay vector on the segment ray changes with position. In this way, the ray transmission process and two-way delay information are integrated and filled into the grid points of the area to be imaged to form a two-way delay matrix.

[0038] Step 3: Fill in the default matrix data.

[0039] Because the ray propagation model uses the antenna position as the origin and emits multiple rays within a certain angle, the azimuth spacing between adjacent rays is large when the rays propagate to deeper locations. Therefore, the two-way delay matrix generated by the ray propagation model is a default matrix, meaning that missing values ​​are irregularly present within the matrix. From a range-azimuth perspective, the two-way delay matrix exhibits characteristics such as dense at the top and sparse at the bottom, dense in the middle and sparse on both sides. Using such a default matrix during imaging results in too few coherent stacks of backscattered echoes from deep targets, resulting in unclear and weak imaging of these targets. To avoid this, the default matrix needs to be padded with data.

[0040] First, the ray propagation range is calibrated according to the two rays with the largest and smallest angles in the ray propagation model. The template matrices with values ​​of 1 and 0 are inside and outside the range respectively. , its significance lies in judging whether an area is a ray propagation area within a certain beam width.

[0041] Then the round-trip delay matrix Perform one-dimensional linear interpolation in the distance direction to fill the two-way delay information at the missing positions in the matrix to obtain the interpolation matrix Since linear interpolation filling is associative, there are linearly related pseudo-delay points outside the ray propagation area, so template matrix matching is required. and The Hadamard product of the matrix is ​​filled :

[0042] ;

[0043] Step 4: Cross-correlation backprojection imaging.

[0044] Because the horizontal distribution of the medium in a layered scene is relatively constant, the ray propagation model corresponding to each measured subaperture is consistent during the measurement process. The two-way delay matrix only needs to be calculated once, significantly reducing computational time. Knowing the two-way delay matrix for each subaperture allows for coherent superposition of backscattered echoes, followed by signal cross-correlation, to perform cross-correlation back-projection (CBP) imaging.

[0045] In the classic BP imaging algorithm, if the number of subapertures is , where the backscatter echo of a single measurement process coherently generates a sub-image, and the process is as follows:

[0046] ；

[0047] in It is sub-images, It is The echo data is obtained by measuring the sub-images. The round-trip delay in When the delay data in corresponds to The sub-image can be obtained by superimposing the amplitude of The sub-images are then superimposed according to the antenna positions to obtain the BP imaging result:

[0048] ；

[0049] For the echo intensity of a single imaging point, its amplitude is determined by Backscatter echo from a sub-aperture However, this simple superposition process does not take into account the correlation between the echo data of each channel, so the cross-correlation operation of the signal is introduced to reduce the influence of side lobes and clutter. The elements are multiplied without repetition and then superimposed. The superposition process can be expressed as:

[0050] ;

[0051] Example

[0052] To verify the proposed ray-theory-based rapid layered backprojection imaging method for ground-penetrating radar (GPR) in this paper, a simulation experiment was designed for a layered medium. Using the electromagnetic simulation software GPRMAX, the simulation employed a Ricker wavelet with a center frequency of 1.5 GHz. The transmitting and receiving antennas were positioned close to the surface, with a stepping distance of 1 cm and a time window of 20 ns. The acquired B-scan data consisted of 101 A-scan channels.

[0053] Scenes such as Figure 2 As shown, in addition to the topmost air layer, there are three underground dielectric layers. The relative dielectric constants of the three layers are 4.0, 8.0 and 6.0 respectively, and the thicknesses are 30 cm, 30 cm and 40 cm respectively. The relative magnetic permeability of each layer is 1, and the conductivity is 0.001 S / m. In the scene, a metal cylindrical target with a radius of 2 cm is placed on the center of each layer. The Ray-theory Layered Cross-correlation Back-Projection (RLCBP) algorithm based on ray theory proposed in the present invention is used for imaging, and two frequency domain algorithms suitable for layered media, Phase Shift Migration (PSM) and Multilayer Range Migration (MRM), are used for comparative verification, and the imaging results are as follows: Figure 3 As shown in the figure, the computation times for the RLCBP, PSM, and MRM algorithms are 1.1168s, 1.3654s, and 1.2659s, respectively. The RLCBP algorithm overcomes the high computational complexity commonly associated with the BP algorithm. Furthermore, the imaging results demonstrate that the signal cross-correlation operation effectively suppresses sidelobes and clutter, a processing capability not available in frequency-domain algorithms such as MRM and PSM. Compared to the other two algorithms, RLCBP offers better focusing performance for deep targets and higher imaging resolution.

[0054] The above specific embodiments merely illustrate the design principles of the present invention. The shapes and names of the components described herein may vary and are not limiting. Therefore, those skilled in the art may modify or substitute equivalents for the technical solutions described in the above embodiments. Such modifications and substitutions, without departing from the inventive spirit and technical solutions of the present invention, shall fall within the scope of protection of the present invention.

Claims

1. A layered cross-correlation backprojection imaging method for ground penetrating radar based on ray theory, characterized in that: The method comprises the following steps: Step 1: Establish a ray propagation model in layered media: Based on Snell's law, establish a ray propagation model in layered media under a common offset system; Step 2: Generate a round-trip delay matrix based on the ray propagation model: The round-trip delay based on the propagation model is calculated for coherent superposition in the subsequent imaging process. Since the round-trip delay of an imaging point in a homogeneous medium is positively correlated with the distance, the round-trip delay information of any point on the transmission path can be quickly calculated based on the ray propagation model and the round-trip delay at the refraction point. Step 3, default matrix data filling: In the entire area to be imaged, only the round-trip delay information on the transmission path is known. It has the characteristics of dense top and sparse bottom, dense in the middle and sparse on both sides. In order to avoid the low energy caused by too few coherent superposition times of deep targets during imaging, the round-trip delay matrix is ​​filled with default matrix; Step 4: Cross-correlation backprojection imaging: The filled time delay matrix provides echo energy information during the imaging point superposition process. Back-projection coherent superposition is performed on the entire imaging area, and signal cross-correlation operation is introduced to suppress sidelobe and clutter energy to obtain the final high-precision imaging result.

2. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 1, characterized in that: In the step 1, The propagation model of the ray is regarded as the signal propagation path of a ground penetrating radar measurement. The antenna position is The signal propagation path It can be expressed as: ； in Representative starting point The antenna transmits Ray trajectory, It is the angle between the ray direction and the antenna measurement direction.

3. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 2, characterized in that: In step 2, solve Round-trip delay at the end point of the segment ray The method is: ； ； in Indicates that electromagnetic waves The propagation speed in the layer medium is such that the electromagnetic wave needs to go back and forth twice on the path to reach each point. Therefore, when calculating the delay between transmission and reception, the wave speed needs to be set to half of the original one. Since the medium of each layer is uniform, the distance of the sub-rays in this layer is directly proportional to the round-trip delay.

4. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 1, characterized in that: In step 2, the step length Make this request: ； in is the wave velocity in the medium, Indicates the minimum step size of the imaging grid points.

5. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 3, characterized in that: In step 2, the physical position of any discrete point on the entire ray propagation path and the corresponding round-trip delay are used to convert the ray propagation model into a round-trip delay matrix of the imaging grid: ； in It is a two-dimensional matrix data, indicating that the coordinates in the imaging grid matrix are Round trip delay information at , Indicates the first The first track The two-way delay vector on the segment ray varies with position.

6. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 1, characterized in that: In step 3, the ray propagation range is calibrated according to the two rays with the largest and smallest angles in the ray propagation model, and the template matrices with values ​​of 1 and 0 are respectively inside and outside the range. , its significance lies in judging whether an area is a ray propagation area within a certain beam width.

7. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 1, characterized in that: In step 3, the round trip delay matrix Perform one-dimensional linear interpolation in the distance direction to fill the two-way delay information at the missing positions in the matrix to obtain the interpolation matrix Since linear interpolation filling is correlated, there are linearly correlated pseudo-delay points outside the ray propagation area, so template matrix matching is required; and The Hadamard product of the matrix is ​​filled : 。 8. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 1, characterized in that: In step 4, since the medium distribution in the horizontal direction of the layered scene is relatively unchanged, the ray propagation model corresponding to each measured sub-aperture is consistent during the measurement process, and the round-trip delay matrix only needs to be calculated once.

9. The ray theory-based ground penetrating radar layered cross-correlation backprojection imaging method according to claim 1, characterized in that: In step 4, the amplitude of the echo intensity of a single imaging point is determined by Backscatter echo from a sub-aperture However, this simple superposition process does not take into account the correlation between the echo data of each channel, so the cross-correlation operation of the signal is introduced to reduce the influence of side lobes and clutter. For a single imaging point, this The elements are multiplied without duplication and then superimposed. From the perspective of sub-images, the superposition process can be expressed as: 。

Images