Ground penetrating radar B-Scan atlas compressed sensing imaging method based on rotation detection

Abstract

The invention discloses an imaging method suitable for a ground penetrating radar (GPR) rotation detection scene. The method comprises the following steps: constructing a sparse matrix strictly conforming to an actual physical model, and performing physical modeling only for a beam range corresponding to a rotation azimuth angle of a beam center of the ground penetrating radar, so that the method has good interpretability; a Gaussian random matrix meeting a finite equidistant criterion is adopted as a projection matrix; taking a reverse Huber norm as a target function, and establishing a convex optimization model of an imaging problem; efficiently solving the optimization problem by using a solver with a primal-dual interior point method; in addition, an evaluation index comprehensively considering an imaging reflection coefficient and a reconstruction signal error is provided, and imaging quality evaluation is carried out by jointly adopting a root mean square error (RMSE) and a Hamming distance. The method can be applied to multi-target simulation and physical imaging scenes based on GPR rotation detection, the imaging quality and the result reliability can be effectively improved, and the development of the ground penetrating radar detection technology is promoted.

Description

Technical Field

[0001] The present invention belongs to the field of digital signal processing, and in particular relates to a ground penetrating radar B-Scan spectrum compressed sensing imaging method based on rotation detection. Background Art

[0002] Among the many existing GPR imaging algorithms, the time domain backprojection algorithm is the most intuitive and simple, but it has the disadvantages of large computational complexity and high sidelobe and clutter energy; the frequency domain range offset algorithm is only effective for imaging scenarios with one layer of medium and fails in actual detection scenarios; the inverse scattering imaging algorithm has too much computational complexity and requires high precision to establish the scene model, which seriously limits its practical application.

[0003] Given the high sampling rate, large data volume and difficulty in real-time processing of ultra-wideband signals, Donoho, Candes, Romberg, Tao and others proposed the compressed sensing (CS) theory in 2006. This is a new signal reconstruction theory that does not rely on the Shannon-Nyquist sampling theorem. It provides ideas for breaking through the many limitations of traditional sampling on ultra-wideband signal processing, and greatly facilitates signal acquisition, storage, transmission and processing.

[0004] CS shows that as long as the signal is compressible or sparse in a certain transform domain, the high-dimensional signal obtained by the transform can be projected onto a low-dimensional space using a measurement matrix that is independent of the sparse basis. Then, by solving an optimization problem, the original signal can be reconstructed with high probability from these small projections. Generally, potential targets such as underground cavities and buried metal objects occupy only a small portion of the entire target space. Therefore, the prior knowledge about GPR images is that they are sparse targets.

[0005] Regarding the application of CS to ground-penetrating radar (GPR) imaging, Gurbuz et al. in the United States first applied compressed sensing theory to GPR imaging and found that, under the same conditions, compressed sensing imaging quality far exceeded that of traditional time-domain and frequency-domain algorithms, achieving high-resolution, low-sidelobe clutter, and robustness to various interference effects. Qu et al. proposed a compression method, but did not consider nonlinearities and pathological conditions, requiring extensive computation during reconstruction. Inspired by CS, to improve data acquisition speed and achieve accurate reconstruction of each A-scan signal, ABSuksmono proposed a compressed SFCW GPR system. They demonstrated the applicability of the proposed system through simulations of single-cycle waveforms and actual measurements using a vector network analyzer within a GPR test range. Summary of the Invention

[0006] Purpose of the Invention: The application of compressed sensing to imaging of common-center and common-offset ground-penetrating radar (GPR) data is becoming increasingly mature. However, this team has proposed a new method for GPR rotation observation for advanced detection. The hardware equipment of this new method's GPR rotation observation system includes a spatially compensated three-dimensional radar antenna, a high-precision three-dimensional radar bracket, a shielded transmission line, a joystick, and a horizontal rotation frame. The radar's pitch angle is adjusted by rotating the joystick, and the horizontal angle is adjusted by turning the horizontal rotation frame. The coordinated horizontal and pitch angles enable the radar to perform comprehensive and intensive detection throughout the entire space. Therefore, based on the basic framework of compressed sensing, the present invention adaptively performs profile imaging on rotating GPR.

[0007] To achieve the above purpose, the present invention adopts the following technologies:

[0008] A ground penetrating radar B-Scan image compression sensing imaging method based on rotation detection includes the following steps:

[0009] (1) Constructing a sparse representation matrix: Based on the main beam pointing direction of the ground penetrating radar antenna at each azimuth angle during its rotation, the corresponding spatial coverage range is selected, and combined with the radar wave propagation path and time delay relationship, a sparse matrix reflecting the physical propagation characteristics of the radar is constructed to achieve physical modeling of the target area and ensure that the model has physical interpretability;

[0010] (2) Generate projection matrix: Use a Gaussian distribution random matrix that satisfies the finite isometric property as the projection matrix;

[0011] (3) Constructing the optimization objective function: The imaging problem is formulated as a convex optimization problem, and the more robust inverse Huber function is introduced as the objective function;

[0012] (4) Solving the optimization problem: Using the primal-dual interior point algorithm to numerically solve the above optimization model, obtain the optimal sparse solution and obtain the imaging result;

[0013] (5) Quality evaluation indicators: Starting from the two dimensions of the accuracy of the imaging reflection coefficient and the similarity of the time domain waveform, RMSE is used to evaluate the pixel value difference of the time domain waveform, and Hamming distance is used to measure the structural difference of the reflection coefficient binary image to comprehensively evaluate the imaging effect.

[0014] Furthermore, the specific steps of constructing the sparse matrix in step 1 include:

[0015] Step 1: Parameter initialization, set the receiving antenna position (rx_pos) and the transmitting antenna position (tx_pos), and initialize the sparse matrix Ψ with a zero matrix;

[0016] Step 2: Angle traversal, iterating each rotation angle of the ground penetrating radar (that is, each angle in the rotation angle range array);

[0017] Step 3: Define the beam range, calculate the antenna position (x0 = rx_pos, y0 = tx_pos), and calculate the beam boundary angle θ with a beam width of 80° left and θ right ;

[0018] Step 4: Grid generation, creating a two-dimensional space grid corresponding to the coordinates in the x and y directions;

[0019] Step 5: Determine the points in the beam. According to the direction angle of the coordinate point relative to the radar position, select the points in the beam angle range [θ left ,θ right ] grid points within;

[0020] Step 6: Delay calculation, for all valid grid points, calculate the round-trip propagation time delay τ = (τ rx +τ tx ) / 2;

[0021] Step 7: Signal interpolation: Use the interp1 function of MATLAB to interpolate the A-scan signal s on the calculated time delay τ;

[0022] Step 8: Update the sparse matrix, normalize the interpolated signal with the l2 norm, and fill its value into the corresponding position of the sparse matrix Ψ.

[0023] Furthermore, the projection matrix Φ defined in step 2 is as follows:

[0024]

[0025] In order to meet the finite equidistance criterion, a Gaussian random matrix is ​​used instead, where N A Represents the total rotation of the ground penetrating radar N A Angles.

[0026] Furthermore, the reverse Huber function defined in step 3 is as follows:

[0027]

[0028] Among them, α and ρ are adjustable parameters, N x 、N y The number of grid points that the detection profile is discretized into at equal intervals of Δx and Δy in the azimuth and range directions, are the reflection coefficients at the corresponding imaging grid points respectively, and function B is the inverse Huber function, which is defined as follows:

[0029]

[0030] The definition of m is as follows:

[0031]

[0032] Where y is the original B-Scan atlas data, assuming there is a set of base So that y is is sparse, that is, there is a Make is k-sparse (i.e. There are at most k non-zero elements).

[0033] Furthermore, the quality evaluation index defined in step 5 is as follows:

[0034]

[0035] Among them, the definitions of isreflect, RMSE and Hamming function are as follows:

[0036]

[0037] Where H*W represents the size of the X and Y maps; numel is a function of the total number of statistical elements.

[0038] Compared with the prior art, the present invention has the following advantages:

[0039] 1. Breaking through the conventional linear detection model, this paper applies compressed sensing technology to rotating ground-penetrating radar imaging for the first time. Existing research has primarily focused on GPR profile imaging based on linear propulsion, failing to effectively adapt to the nonlinear coverage data collected by rotating antennas. This invention exploits the characteristics of antenna rotation detection by constructing a sparse representation matrix that dynamically matches the beam direction and propagation path, enabling physical modeling of rotating detection data, significantly improving the applicability and accuracy of imaging methods.

[0040] 2. An innovative dual-index system for imaging quality is proposed, integrating RMSE and Hamming distance. This dual-index system considers both the pixel-level error of the time domain waveform (RMSE) and the structural similarity of the reflectance coefficient map (Hamming distance), filling a gap in existing GPR compressed sensing imaging quality evaluation standards.

[0041] 3. A robust optimization objective model based on the inverse Huber function was constructed, effectively improving the stability and convergence of the reconstruction algorithm in environments with strong noise and abnormal interference. The inverse Huber function combines the advantages of both the L1 and L2 norms, offering stronger interference resistance. Compared to traditional L1 norm optimization models, it offers superior reconstruction quality and convergence speed, making it more suitable for GPR multi-target imaging tasks in complex underground scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 This is a flow chart of the compressed sensing imaging method of the ground penetrating radar B-Scan spectrum based on rotation detection of the present invention;

[0043] Figure 2 A rotation detection scene diagram of the ground penetrating radar of the present invention;

[0044] Figure 3 Ricker wavelet diagrams of different center frequencies of the present invention;

[0045] Figure 4 (a) Original reflection coefficient distribution diagram of the present invention; (b) Original B-Scan signal;

[0046] Figure 5 The detection area of ​​the ground penetrating radar of the present invention at different angles, with the beam width set to 60° (the black grid is the model space area, and the black circle grid is the detection area) (a) The antenna rotation angle is at 0°; (b) The antenna rotation angle is at 90°; (c) The antenna rotation angle is at 180°; (d) The antenna rotation angle is at 270°;

[0047] Figure 6 This is a relationship diagram of the solution matrix of the compressed sensing of the ground penetrating radar of the present invention;

[0048] Figure 7 The imaging results of different algorithms in a simulation example of the present invention are as follows: (a) LASSO; (b) ADMM; (c) ISTA; (d) Ours. DETAILED DESCRIPTION

[0049] The technical solution of the present invention is further described below with reference to the accompanying drawings, but is not limited thereto. Any modification or equivalent replacement of the technical solution of the present invention that does not depart from the spirit and scope of the technical solution of the present invention should be included in the scope of protection of the present invention.

[0050] The specific operation process of the ground penetrating radar B-Scan atlas compressed sensing imaging method based on rotation detection is as follows: Figure 1 As shown, the following steps are included:

[0051] (1) Construct a data model, and use pulsed ground penetrating radar to perform rotational detection of the entire space. The detection scene is as follows: Figure 2 As shown. Assume that the antenna rotates horizontally at equal intervals of Δθ. Let the ground penetrating radar transmit pulse be s(t). Considering the time delay and waveform expansion of the target echo model at the detection profile grid point q of the ground penetrating radar at the pth angle in the horizontal direction, the receiving antenna signal model can be established as:

[0052] s p (t) = r q s(t-τ(p,q))

[0053] The corresponding data is stored in matrix form:

[0054]

[0055] Where τ(p,q) is the two-way delay of the electromagnetic wave from the pth angle to the midpoint q of the detection area profile, r q is the scattering coefficient at point q; F s is the time domain sampling frequency, N s is the number of sampling points of s(t) in the time domain.

[0056] The ground penetrating radar rotates along the horizontal direction for detection, with a total rotation of N A N A The two-dimensional B-Scan echo data matrix composed of A-Scan is recorded as s p The column vector is an A-Scan data of the ground penetrating radar at an angle p in the horizontal direction. Arrange each data in S in a column and record it as

[0057] When the ground penetrating radar rotates to the pth angle, the two-way delay τ(p,q) from the grid point q is A ], the position of point p in the detection profile grid is recorded as row i, column j, i∈[1,N y ],j∈[1,N x ]; N x 、N y The number of grid points into which the detection profile is discretized at equal intervals of Δx and Δz in the azimuth and distance directions respectively; the propagation speed of electromagnetic waves in a uniform medium is taken as c is the speed of light, ε r is the dielectric constant of the propagation medium.

[0058] Here, the pulse signal emitted by the GPR selected in the present invention is the Ricker wavelet, such as Figure 3 shown.

[0059] If feasible, design relevant simulation models. Figure 4As shown in (a), a two-dimensional detection area with a length of 15m and a width of 10m is established. Two rectangular target cross-section holes are set in the area. The relative dielectric constant of the background medium is 2.5. The antenna is located at (5m, 12m). Figure 4 (b) shows the raw B-scan data obtained from the rotational detection (only the mean subtraction operation was performed to remove direct waves). The detection method is to rotate the antenna from 0° to 180°, completing the detection in 1° intervals. A total of 181 A-scans are collected. The spatial grid spacing is 0.1m × 0.1m, and the model area has a total of 150 × 100 grid points. The transmitter is a Ricker wavelet with a main frequency of 100MHz.

[0060] (2) Construct the sparse matrix Ψ. The specific construction steps are as follows:

[0061] Step 1: Based on the B-Scan data obtained above, perform parameter initialization, set the positions of the transmitting antenna (tx_pos) and the receiving antenna (rx_pos), and initialize the sparse matrix Ψ with a zero matrix;

[0062] Step 2 Angle traversal, traverse each angle in the rotation angle array in turn, simulating the radar observation at different angles;

[0063] Step 3: Define the beam range. Calculate the corresponding beam boundary angle θ based on the antenna position (x0 = rx_pos, y0 = tx_pos) and the set 80° beam width. left and θ right ;

[0064] Step 4: Grid generation, generating a two-dimensional space coordinate grid covering the x-axis and y-axis;

[0065] Step 5: Determine the points within the beam. According to the azimuth angle of each grid point relative to the radar position, select the points that fall within the beam range [θ left ,θ right ] within the valid points;

[0066] Step 6 Delay calculation: For all valid grid points, calculate the round-trip propagation time delay τ = (τ rx +τ tx ) / 2;

[0067] Step 7: Signal interpolation: Use the interp1 function in MATLAB to interpolate the A-scan signal s based on the calculated time delay τ.

[0068] Step 8: Update the sparse matrix. Normalize the interpolation result with the l2 norm and fill it into the corresponding position of the sparse matrix Ψ.

[0069] Among them, according to a certain beam width, the sparse matrix Ψ corresponding to the ground penetrating radar at the pth rotation angle is established p for:

[0070]

[0071] like Figure 5 As shown in the figure, only grid points within the beam width will be assigned values, and those outside the range will all be set to 0.

[0072] An A-Scan echo s when the ground penetrating radar rotates to the pth angle p The sparse matrix Ψ can be used p Expressed as:

[0073]

[0074] in, are the reflection coefficients at the corresponding imaging grid points. p =Φ1s p =Φ1Ψ p r, where Φ1 is the observation matrix.

[0075] Globally, the data collected by the compressed sensing end corresponds to the low-dimensional Y obtained by linear projection of the B-Scan through the projection matrix Φ, which is mathematically expressed as:

[0076]

[0077] in,

[0078]

[0079] Group the column vectors The column vectors are arranged in a column to form a large column vector, which is recorded as

[0080] According to the compressed sensing theory, ,A cs The requirement is to meet the finite equidistance criterion, and the Gaussian random matrix Φ can better meet the requirement. The reflection coefficient r at each grid point can be solved by the norm minimization problem, that is,

[0081]

[0082] The two-dimensional imaging results can be expressed as a two-dimensional matrix. Each element in the matrix corresponds to the energy at the imaging grid. Converting the solved column vector r into a two-dimensional matrix is ​​the imaging result.

[0083] In summary, it can be concluded that Figure 6The matrix relationship of compressed sensing in the solution process is shown.

[0084] (3) Construct the optimization objective function. For the original B-Scan atlas y to be imaged, assume that there is a set of bases So that y (column vector) is is sparse, ie st is k-sparse (i.e. There are at most k non-zero elements), the expression corresponding to y is as follows:

[0085]

[0086] Generally speaking, the present invention "measures" not the real image y, but its transformed form, denoted as m. The present invention assumes a linear measurement model, using the matrix Its mathematical expression is:

[0087]

[0088] The objective function optimized by the present invention is based on the reverse Huber norm, as follows:

[0089]

[0090] Among them, α and ρ are adjustable parameters, and function B is the inverse Huber function, which is defined as follows:

[0091]

[0092] (4) In order to solve the above optimization problem, the present invention uses the SDPT3 solver (i.e., the primal-dual interior point method) built into the CVX toolbox to solve and optimize the function.

[0093] (5) In order to test the imaging effect, the present invention proposes a joint error index based on the solved sparse coefficient r and the reconstructed signal (Ψr), involving RMSE and Hamming distance:

[0094]

[0095] Among them, the definitions of isreflect, RMSE and Hamming function are as follows:

[0096]

[0097] Where H*W represents the size of the X and Y maps; numel is a function of the total number of statistical elements.

[0098] The imaging indicators of different algorithms on the simulation data are shown in Table 1. The imaging comparison effect is as follows Figure 7 shown.

[0099] Table 1 Imaging index results of different algorithms

[0100]

[0101] The above embodiments further illustrate the purpose, technical solutions and advantages of the present invention in detail. It should be understood that the above embodiments are only preferred implementation methods of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made to the present invention within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A ground penetrating radar B-Scan image compression sensing imaging method based on rotation detection, characterized by The ground penetrating radar rotation detection imaging method is carried out according to the following steps: (1) Constructing a sparse representation matrix: Based on the main beam pointing direction of the ground penetrating radar antenna at each azimuth angle during its rotation, the corresponding spatial coverage range is selected, and combined with the radar wave propagation path and time delay relationship, a sparse matrix reflecting the physical propagation characteristics of the radar is constructed to achieve physical modeling of the target area and ensure that the model has physical interpretability; (2) Generate projection matrix: Use a Gaussian distribution random matrix that satisfies the finite isometric property as the projection matrix; (3) Constructing the optimization objective function: The imaging problem is formulated as a convex optimization problem, and the more robust inverse Huber function is introduced as the objective function; (4) Solving the optimization problem: Using the primal-dual interior point algorithm to numerically solve the above optimization model, obtain the optimal sparse solution and obtain the imaging result; (5) Quality evaluation indicators: Starting from the two dimensions of the accuracy of the imaging reflection coefficient and the similarity of the time domain waveform, RMSE is used to evaluate the pixel value difference of the time domain waveform, and Hamming distance is used to measure the structural difference of the reflection coefficient binary image to comprehensively evaluate the imaging effect.

2. The method for compressed sensing imaging of ground penetrating radar B-Scan images based on rotation detection according to claim 1 is characterized in that The specific steps of constructing the sparse matrix in step 1 include: Step 1: Parameter initialization, set the receiving antenna position (rx_pos) and the transmitting antenna position (tx_pos), and initialize the sparse matrix Ψ with a zero matrix; Step 2: Angle traversal, iterating each rotation angle of the ground penetrating radar (that is, each angle in the rotation angle range array); Step 3: Define the beam range, calculate the antenna position (x0 = rx_pos, y0 = tx_pos), and calculate the beam boundary angle θ with a beam width of 80° left and θ right ; Step 4: Grid generation, creating a two-dimensional space grid corresponding to the coordinates in the x and y directions; Step 5: Determine the points in the beam. According to the direction angle of the coordinate point relative to the radar position, select the points in the beam angle range [θ left ,θ right ] grid points within; Step 6: Delay calculation, for all valid grid points, calculate the round-trip propagation time delay τ = (τ rx +τ tx ) / 2; Step 7: Signal interpolation: Use the interp1 function of MATLAB to interpolate the A-scan signal s on the calculated time delay τ; Step 8: Update the sparse matrix, normalize the interpolated signal with the l2 norm, and fill its value into the corresponding position of the sparse matrix Ψ.

3. The method for compressed sensing imaging of ground penetrating radar B-Scan images based on rotation detection according to claim 1 is characterized in that The projection matrix Φ defined in step 2 is as follows: In order to meet the finite equidistance criterion, a Gaussian random matrix is ​​used instead, where N A Represents the total rotation of the ground penetrating radar N A Angles.

4. The method for compressed sensing imaging of ground penetrating radar B-Scan images based on rotation detection according to claim 1 is characterized in that The reverse Huber function defined in step 3 is as follows: Among them, α and ρ are adjustable parameters, N x 、N y The number of grid points that the detection profile is discretized into at equal intervals of Δx and Δy in the azimuth and range directions, are the reflection coefficients at the corresponding imaging grid points respectively, and function B is the inverse Huber function, which is defined as follows: The definition of m is as follows: Where y is the original B-Scan atlas data, assuming there is a set of base So that y is is sparse, that is, there is a Make is k-sparse (i.e. There are at most k non-zero elements).

5. The method for compressed sensing imaging of ground penetrating radar B-Scan images based on rotation detection according to claim 1 is characterized in that The quality evaluation index defined in step 5 is as follows: Among them, the definitions of isreflect, RMSE and Hamming function are as follows: Where H*W represents the size of the X and Y maps; numel is a function of the total number of statistical elements.

Images