Odometer estimation method based on ground penetrating radar

Abstract

The invention discloses a speedometer estimation method of a ground penetrating radar. The method comprises the following steps: carrying out background removal, DC offset removal, triangular band-pass filtering, SEC gain compensation, wavelet denoising and Gaussian smoothing processing on GPR original signals in a public data set to obtain multiple frames of B scanning images; an original odometer estimation model is constructed and trained, and a trained odometer estimation model is obtained; and acquiring GPR original signals of a series of underground ground features by using a ground penetrating radar mounted on a traveling vehicle or robot, processing the GPR original signals into B scanning images, inputting two adjacent frames of B scanning images into the trained speedometer estimation model, and outputting the traveling distance of the vehicle or robot. According to the speedometer estimation model disclosed by the invention, the estimation precision of the Euclidean distance between two adjacent frames of B scanning images can be remarkably improved by jointly extracting the difference characteristics and the similarity characteristics of the B scanning images, and the speedometer estimation model is particularly suitable for severe weather or scenes with stable underground structures, such as tunnels, parking lots and the like.

Description

Technical Field

[0001] The present invention belongs to the technical field of ground penetrating radar scanning image detection, and in particular relates to an odometer estimation method based on ground penetrating radar. Background Art

[0002] When using sensors to measure the distance traveled by a vehicle or robot, the most commonly used is the wheel encoder. However, it has a drawback: if the wheel slips or gets stuck, the measurement value of the wheel encoder will be invalid.

[0003] Accurately positioning robots and vehicles in challenging weather and environmental scenarios is crucial for autonomous driving. Currently, commonly used positioning methods rely on onboard sensors such as global positioning system (GPS) receivers, cameras, lidar, and inertial measurement units (IMUs), mainly including the following methods:

[0004] ① Learning-based GPR odometry model: Using a residual network (ResNet-18), feature maps are extracted from two B-scan images with significant features. A sliding window is then used to compare the similarity of the two sets of feature maps to obtain a set of most similar window indices. The window indices are then put into a fully connected layer for regression to obtain an estimate of the distance traveled by the vehicle or robot.

[0005] ②CNN-LSTM model: Multiple pairs of consecutive adjacent B-scan images are placed into a convolutional neural network (CNN) to obtain a set of feature vectors, and then this set of feature vectors is placed into a long short-term memory network (LSTM) to obtain an estimated value of the distance traveled by the vehicle or robot.

[0006] ③ Multimodal odometry model: By extracting key reflection peak points from the B-scan image, a subsurface feature matrix (SFM) that can characterize the spatial distribution of underground features is constructed. Then, the optimal horizontal matching offset between the SFMs at adjacent moments is calculated and converted into physical distance to obtain an estimate of the distance traveled by the vehicle or robot.

[0007] However, these methods are severely hampered in environments such as urban canyons and tunnels, or in adverse weather conditions, posing significant safety challenges for autonomous driving. Furthermore, they rely on artificial feature design, which is highly subjective and computationally complex, and fail to detect highly similar adjacent B-scan images.

[0008] Ground-penetrating radar (GPR) is a less-utilized sensor that effectively detects underground features. Compared to wheel encoders, GPR data acquisition is unaffected by wheel slip. In contrast, the underground structure and features of urban roads remain stable and are less affected by these adverse conditions, providing new opportunities for improving the reliability and robustness of positioning in adverse weather or environmental conditions. However, existing GPR odometry technology suffers from large distance estimation errors when processing adjacent B-scan images due to the minimal feature differences between them. Summary of the Invention

[0009] In order to solve the above problems, the present invention aims to provide an odometer estimation method based on ground penetrating radar.

[0010] To achieve the above-mentioned object, the method for estimating odometer based on ground penetrating radar provided by the present invention comprises the following steps performed in sequence:

[0011] 1) Using the publicly available CMU-GPR dataset, the GPR raw signals were subjected to background removal, DC offset removal, triangular bandpass filtering, SEC gain compensation, wavelet denoising, and Gaussian smoothing to obtain multi-frame B-scan images;

[0012] 2) constructing an original odometry estimation model, and then inputting any two adjacent B-scan images obtained in step 1) into the original odometry estimation model for training to obtain a trained odometry estimation model;

[0013] 3) Using a ground-penetrating radar (GPR) mounted on a moving vehicle or robot, a series of GPR raw signals of underground features are collected. These signals are then processed according to the method in step 1) to obtain multiple B-scan images. Two adjacent B-scan images are then fed into the trained odometry estimation model to output the distance traveled by the vehicle or robot.

[0014] In step 1), the method of performing background removal, DC offset removal, triangular bandpass filtering, SEC gain compensation, wavelet denoising and Gaussian smoothing on the GPR original signal to obtain a multi-frame B-scan image is:

[0015] 1.1) Background removal

[0016] Calculate the average value of all scan signals on the horizontal track and subtract it from the GPR original signal to obtain the background removal component. The formula is as follows:

[0017]

[0018] Where N is the number of horizontal channels, t is the time sampling point, x is the horizontal channel number, and s(t, x) is the xth row of the original signal;

[0019] 1.2) DC offset removal

[0020] Perform cubic polynomial curve fitting on each scan trace to obtain the baseline component:

[0021]

[0022] The baseline component is then subtracted from the GPR original signal s(t) to obtain the component after removing the DC offset, as shown in the following formula:

[0023] s dewow (t)=s(t)-s base (t)

[0024] Among them, c k are the polynomial coefficients;

[0025] 1.3) Triangular bandpass filtering

[0026] The filter coefficients are determined by the following equation:

[0027]

[0028] Among them, f min =200MHz, f max =850MHz is the passband boundary, f s is the sampling frequency, w tri is a triangular window function, n is an integer variable, and its value range is -K≤n≤K, where K=(N-1) / 2, and N is the filter length or order;

[0029] Based on the filter coefficients and the GPR original signal s(t), the filtered components are calculated as follows:

[0030] s filtered (t)=s(t)*h[n]

[0031] 1.4)SEC gain compensation processing

[0032] First, define the time-varying gain function:

[0033] G(t)=t b ·e a·t

[0034] Among them, a is the exponential gain coefficient; b is the power law gain coefficient; then the above time-varying gain function is used to calculate the value of the time threshold value t th The GPR original signal s(t) is gain saturated to obtain the time-varying gain function G after limiting clip (t):

[0035]

[0036] Based on the time-varying gain function G after limiting clip (t) and the GPR original signal s(t), and the enhanced component is calculated:

[0037] s sec (t)=s(t)·G clip (t)

[0038] 1.5) Wavelet denoising

[0039] The GPR original signal s(t) is decomposed based on the Daubechies wavelet basis function, and the formula is as follows:

[0040] c j,k = <s,ψ j,k >

[0041] The Bayesian adaptive threshold algorithm is used to determine the optimal threshold:

[0042]

[0043] in, is the noise variance estimate;

[0044] Based on the Daubechies wavelet basis function and the optimal threshold, the wavelet denoising component is obtained:

[0045]

[0046] 1.6) Gaussian smoothing

[0047] The Gaussian kernel function with a standard deviation of σ is used to perform spatial smoothing on the GPR original signal s(t). The formula is as follows:

[0048]

[0049] The time domain smoothing is achieved through convolution operation to obtain the smoothed component:

[0050]

[0051] In step 2), the original odometer estimation model is constructed, and then any two adjacent frames of B-scan images obtained in step 1) are input into the original odometer estimation model for training. The method for obtaining the trained odometer estimation model is:

[0052] 2.1) Constructing a raw odometry estimation model; the raw odometry estimation model consists of a cascade of a feature extraction module, a parallel difference detection module, a similarity detection module, and a fully connected regression module; wherein the feature extraction module consists of a cascade of a backbone network and an adaptive average pooling layer; the backbone network uses the first six layers of the ResNet-18 network from the initial convolutional layer to the fourth residual block; the difference detection module consists of a cascade of a tensor subtraction layer, a 3×3 convolutional layer, a batch normalization layer, an activation function, and a dropout unit; the similarity detection module consists of a cascade of a flattening layer and a cosine similarity calculation unit; and the fully connected regression module consists of two parallel fully connected layers;

[0053] 2.2) Inputting any two adjacent B-scan images obtained in step 1) into the feature extraction module in the above-mentioned original odometry estimation model to extract features and obtain global features;

[0054] First, the backbone network in the feature extraction module is used to extract the B-scan images B at time t-1 and time t through stacking operations of convolution, batch normalization, and ReLU activation. t-1 and B t The spatial features in the 1-channel 224×224 grayscale B-scan image are converted into a 512-channel 7×7 spatial feature map. The formula is as follows:

[0055] F s =ReLU(BN(Conv 3×3 (I)))+Identity(I)

[0056] Where I∈R B×1×224×224 is the batch of B-scan images, ReLU is the ReLU activation function, Identity(I) is the residual connection, Conv is the convolution operation, and BN is the batch normalization operation;

[0057] Then, the adaptive average pooling layer is used to globally compress the above spatial feature map. The formula is as follows:

[0058]

[0059] Among them, H is the height of the spatial feature map, W is the width of the spatial feature map, and c is the channel index;

[0060] Get the global feature F g ∈R B×512×1×1 ;

[0061] 2.3) Using a difference detection module to perform difference detection on the global features obtained in step 2.2) above to obtain difference features;

[0062] The tensor subtraction layer in the difference detection module is used to calculate two global features F g,t-1 and F g,t The absolute difference between them is used to obtain the difference vector representing the feature difference. The formula is as follows:

[0063] ΔF=|F g,t -F g,t-1 |

[0064] Based on the above difference vector ΔF, a 3×3 convolutional layer is used to extract the local difference features of the two global features. The formula is as follows:

[0065]

[0066] Where W∈R 512×512×3×3 is a learnable convolution kernel, represents the convolution operation, b c is the bias term; then based on the above local difference features, the batch normalization layer is used to stabilize the training process and obtain the normalized features. The formula is as follows:

[0067]

[0068] Among them, μ B ,σ B is the batch statistic, ε is a very small constant to prevent division by zero, and β is a learnable bias parameter;

[0069] Finally, based on the above normalized features, we use the activation function and the discard unit to introduce nonlinear activation and regularization to obtain the difference features. The formula is as follows:

[0070] D act =max(0.2×D norm , D norm )

[0071] D t-1,t =Dropout 0.2 (D act )

[0072] Among them, the 20% dropout rate Dropout improves the generalization of the model through random inactivation;

[0073] 2.4) Using a similarity detection module to perform similarity detection on the spatial feature map obtained in step 2.2) above to obtain similarity features;

[0074] The spatial feature map F obtained in step 2.2) is transformed into s Flattened, the formula is as follows:

[0075] φ:R B×512×7×7 →RB×512×49

[0076] Then, the cosine similarity calculation unit is used to calculate the spatial feature map F channel by channel. s The cosine similarity of , continues to obtain the similarity feature, the formula is as follows:

[0077]

[0078] 2.5) Using a fully connected regression module, regress and fuse the difference features obtained in step 2.3) and the similarity features obtained in step 2.4) to obtain a trained odometry estimation model;

[0079] First, the above differential features D t-1，t and similarity feature S t-1，t Independent regression, the formula is as follows:

[0080]

[0081] Among them, W d ,W s ∈R 512 is the learnable weight, b d , b s is the learnable bias value; is used to quantify the absolute intensity of feature changes, It is used to measure the consistency of feature direction; finally, the above absolute strength is adjusted by adaptive weight α and consistency Perform dynamic fusion to obtain the travel distance O of the vehicle or robot t-1，t , the formula is as follows:

[0082]

[0083] Finally, the training is completed to obtain the trained odometry estimation model.

[0084] The odometer estimation method based on ground penetrating radar provided by the present invention has the following beneficial effects:

[0085] The odometry estimation model used can significantly improve the estimation accuracy of the Euclidean distance between two adjacent B-scan images by jointly extracting the difference and similarity features of B-scan images. It is particularly suitable for scenarios with severe weather or stable underground structures, such as tunnels and parking lots. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1 This is a diagram of the odometry estimation model (GPR-OdomNet) in the present invention.

[0087] Figure 2(a) and (b) are two adjacent B-scan images in the present invention; Figure 2 (c) and (d) show the global features F extracted by the feature extraction module of the B-scan image at time t-1 and time t, respectively. t-1 and F t ; Figure 2 (e) and (f) are similarity features and difference features, respectively.

[0088] Figure 3 Schematic diagram of the factor graph model in the present invention.

[0089] Figure 4 This is the comparison result of the odometer performance in the sensor fusion of the present invention.

[0090] Figure 5 The root mean square error (RMSE) analysis diagram of the absolute trajectory error (ATE) over time for different data sets. DETAILED DESCRIPTION

[0091] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0092] The method for estimating odometer based on ground penetrating radar provided by the present invention comprises the following steps performed in sequence:

[0093] 1) Using the publicly available CMU-GPR dataset, the GPR raw signals were subjected to background removal, DC offset removal, triangular bandpass filtering, SEC gain compensation, wavelet denoising, and Gaussian smoothing to obtain multi-frame B-scan images;

[0094] 1.1) Background removal

[0095] Due to antenna coupling and system fixed pattern noise in the GPR raw signal, there is horizontal stripe interference in the image. To eliminate this noise and improve image uniformity, the GPR raw signal s(t) needs to be background-subtracted. This is done by calculating the average value of all scan signals on the horizontal track and subtracting this average value from the GPR raw signal to obtain the background-subtracted component. The formula is as follows:

[0096]

[0097] Where N is the number of horizontal channels, t is the time sampling point, x is the horizontal channel number, and s(t, x) is the xth row of the original signal;

[0098] 1.2) DC offset removal

[0099] In view of the common low-frequency baseline drift problem in GPR systems, this paper uses polynomial fitting technology to correct the GPR raw signal s(t) to effectively remove the DC offset component and retain the true high-frequency signal characteristics. The method is to perform a cubic polynomial curve fitting on each scan track to obtain the baseline component:

[0100]

[0101] The baseline component is then subtracted from the GPR original signal s(t) to obtain the component after removing the DC offset, as shown in the following formula:

[0102] s dewow (t)=s(t)-s base (t)

[0103] Among them, c k are the polynomial coefficients;

[0104] 1.3) Triangular bandpass filtering

[0105] To address the issues of environmental noise interference and signal band distortion, the present invention uses a triangular window bandpass filter to retain the 200-850MHz effective frequency band in the GPR original signal s(t), effectively suppressing out-of-band noise and maintaining the target reflection characteristic frequency band. The filter coefficient is determined by the following formula:

[0106]

[0107] Among them, f min =200MHz, f max =850MHz is the passband boundary, f s is the sampling frequency, w tri is a triangular window function, n is an integer variable, representing the discrete time index position of the filter coefficient, with a value range of -K≤n≤K, where K=(N-1) / 2, N is the filter length or order, used to calculate the filter coefficient value corresponding to each position;

[0108] Based on the filter coefficients and the GPR original signal s(t), the filtered components are calculated as follows:

[0109] s filtered (t)=s(t)*h[n]

[0110] 1.4)SEC gain compensation processing

[0111] In view of the attenuation effect of electromagnetic waves when propagating in the medium, the present invention adopts a composite gain compensation strategy to perform SEC gain compensation on the GPR original signal s(t), thereby significantly enhancing the detectability of deep weak signals. First, define the time-varying gain function:

[0112] G(t)=t b ·e a·t

[0113] Among them, a is the exponential gain coefficient, which is used to compensate for medium absorption; b is the power law gain coefficient, which is used to compensate for spherical diffusion. Then, the above time-varying gain function is used to compensate for the time threshold t th The GPR original signal s(t) is gain saturated to obtain the time-varying gain function G after limiting clip (t):

[0114]

[0115] Based on the time-varying gain function G after limiting clip (t) and the GPR original signal s(t), and the enhanced component is calculated:

[0116] s sec (t)=s(t)·G clip (t)

[0117] 1.5) Wavelet denoising

[0118] To suppress random noise and retain effective signal features, the present invention uses multi-scale wavelet analysis technology to perform wavelet denoising on the GPR original signal s(t), which can suppress noise while maintaining signal edge features. The GPR original signal s(t) is decomposed based on the Daubechies wavelet basis function, and the formula is as follows:

[0119] c j，k =〈s,ψ j,k >

[0120] The Bayesian adaptive threshold algorithm is used to determine the optimal threshold:

[0121]

[0122] in, is the noise variance estimate;

[0123] Based on the Daubechies wavelet basis function and the optimal threshold, the wavelet denoising component is obtained:

[0124]

[0125] 1.6) Gaussian smoothing

[0126] In order to optimize image quality, effectively suppress random fluctuations, and improve the visual quality of the image, the present invention uses a Gaussian kernel function with a standard deviation of σ to perform spatial smoothing on the GPR original signal s(t). The formula is as follows:

[0127]

[0128] The time domain smoothing is achieved through convolution operation to obtain the smoothed component:

[0129]

[0130] Through this multi-stage signal processing, the raw GPR signal s(t) is converted into a high-quality B-scan image, significantly improving the signal-to-noise ratio by 30-50 decibels (primarily due to background removal and wavelet denoising), enhancing deep-layer target detection capabilities by 8-15 times (primarily due to SEC gain compensation processing), ensuring a spatial resolution of 0.015 meters (achieved by triangular bandpass filtering), and optimizing the dynamic range to ±(2-3) times the standard deviation (accomplished through Gaussian smoothing), effectively resolving key technical issues such as noise interference, signal attenuation, and baseline drift in the raw GPR signal s(t). Each processing step is based on a precise mathematical model, and the processing results are highly repeatable and accurate, providing a reliable image foundation for subsurface target identification and interpretation.

[0131] 2) Build an original odometry estimation model, and then input any two adjacent B-scan images obtained in step 1) into the original odometry estimation model for training to obtain a trained odometry estimation model (GPR-OdomNet);

[0132] 2.1) Build Figure 1 The original odometry estimation model shown in FIG4 is composed of a cascade of a feature extraction module, a parallel difference detection module, a similarity detection module, and a fully connected regression module. The feature extraction module is composed of a cascade of a backbone network and an adaptive average pooling layer. The backbone network adopts the first six layers from the initial convolutional layer to the fourth residual block in the ResNet-18 network. The difference detection module is composed of a cascade of a tensor subtraction layer, a 3×3 convolutional layer, a batch normalization layer, an activation function, and a dropout unit. The similarity detection module is composed of a cascade of a flattening layer and a cosine similarity calculation unit. The fully connected regression module is composed of two parallel fully connected layers.

[0133] 2.2) Inputting any two adjacent B-scan images obtained in step 1) into the feature extraction module in the above-mentioned original odometry estimation model to extract features and obtain global features;

[0134] like Figure 2 As shown in (a) and (b), noise from B-scan images is inevitable, which may mask the critical reflection signals. A key step is to suppress noise and non-target reflections to improve the visibility of valid reflection signals.

[0135] First, the backbone network in the feature extraction module is used to extract the B-scan images B at time t-1 and time t through stacking operations of convolution, batch normalization, and ReLU activation. t-1 and B t The spatial features in the 1-channel 224×224 grayscale B-scan image are converted into a 512-channel 7×7 spatial feature map. The formula is as follows:

[0136] F s =ReLU(BN(Conv 3×3 (I)))+Identity(I)

[0137] Where I∈R B×1×224×224 is the batch of B-scan images, ReLU is the ReLU activation function, Identity(I) is the residual connection, Conv is the convolution operation, and BN is the batch normalization operation. This formula retains the original feature information through skip connections and solves the gradient disappearance problem of deep networks. The output spatial feature map F s ∈R B×512×7×7 It contains high-level semantic information and preserves spatial structure.

[0138] Then, the adaptive average pooling layer is used to globally compress the above spatial feature map. The formula is as follows:

[0139]

[0140] Among them, H is the height of the spatial feature map, W is the width of the spatial feature map, and c is the channel index;

[0141] This operation compresses the 7×7 spatial feature map of each channel into a 1×1 scalar to obtain the global feature F g ∈R B ×512×1×1 Its physical significance lies in aggregating spatial information to generate image-level descriptors to eliminate position sensitivity and provide high-level semantic representation for downstream tasks. Figure 2 (c) and (d) show the B-scan images B at time t-1 and time t, respectively. t-1 and B t After the above operation, the global feature F is extracted g,t-1 and F g,t Since the dimension of the extracted global feature is 512, it is not convenient to illustrate. Figure 2 (c) and (d) show global features extracted from one dimension out of every ten dimensions. This step facilitates the effective differentiation and extraction of deep features, including shape and texture information. This feature extraction not only mitigates the interference of noise but also enhances the characteristic representation of the target structure in the image, thereby improving the accuracy of subsequent image analysis.

[0142] 2.3) Using a difference detection module to perform difference detection on the global features obtained in step 2.2) above to obtain difference features;

[0143] Scan image B from B respectively t-1 and B t Extract the global feature F g,t-1 and F g,t After that, the next step is to use these two global features to estimate the travel distance of the vehicle or robot, which requires focusing on the similarities and differences between the two global features. Similarity detection can identify the common points between two global features. However, two adjacent B-scan images are usually very similar, such as Figure 2 (c) and (d). Therefore, it is important to identify the differences between them because they contain important information about the GPR motion.

[0144] The tensor subtraction layer in the difference detection module is used to calculate two global features F g,t-1 and F g,t The absolute difference between them is used to obtain the difference vector representing the feature difference. The formula is as follows:

[0145] ΔF=|F g，t -F g,t-1 |

[0146] Each element in the difference vector ΔF reflects the magnitude of the difference between the two global features along the corresponding feature dimension. Taking the absolute value is to ignore the direction of the difference and only focus on the magnitude of the difference.

[0147] In order to further extract local patterns and structural information, a convolution operation can be used to capture the local correlation within the difference vector ΔF, see Figure 2 (e), for example, identifying whether the differences along certain feature dimensions show spatial regularity. Based on the above difference vector ΔF, a 3×3 convolutional layer is used to extract the local difference features of the two global features, as follows:

[0148]

[0149] Where W∈R 512×512×3×3 is a learnable convolution kernel, represents the convolution operation, b c is the bias term. The key roles of this convolutional layer are: 1) learning the differential correlation between channels; 2) enhancing the response of local variation areas; and 3) maintaining the spatial consistency of features (even though the input is 1×1, the convolution still acts on the channel dimension). Based on the above local differential features, the batch normalization layer stabilizes the training process and obtains the normalized features. The formula is as follows:

[0150]

[0151] Among them, μ B ,σ B is the batch statistic, ε is a very small constant to prevent division by zero, and β is a learnable offset parameter; this operation has three meanings: 1) eliminates internal covariate shift; 2) allows higher learning rates; 3) reduces sensitivity to initialization.

[0152] Finally, based on the above normalized features, we use the activation function and the discard unit to introduce nonlinear activation and regularization to obtain the difference features. The formula is as follows:

[0153] D act =max(0.2×D norm , D norm )

[0154] D t-1,t =Dropout 0.2 (D act )

[0155] Among them, the 20% dropout rate is used to improve the generalization of the model through random inactivation. The output difference feature D t-1,t ∈R B×512×1×1 It encodes the essential characteristics of scene changes, enabling it to capture more complex feature difference patterns.

[0156] 2.4) Using a similarity detection module to perform similarity detection on the spatial feature map obtained in step 2.2) above to obtain similarity features;

[0157] In parallel with difference detection, by calculating the cosine similarity of two spatial feature maps, the similarity features between the two feature maps can be obtained to improve the accuracy of distance estimation, see Figure 2 (f) Directly calculating channel-level cosine similarity on spatial feature maps avoids introducing trainable parameters.

[0158] The spatial feature map F obtained in step 2.2) is transformed into s Flattened, the formula is as follows:

[0159] φ:R B×512×7×7 →R B×512×49

[0160] Then, the cosine similarity calculation unit is used to calculate the spatial feature map F channel by channel. s The cosine similarity of , continues to obtain the similarity feature, the formula is as follows:

[0161]

[0162] The geometric meaning of the above formula is to calculate two spatial feature maps F s，t (c), F s，t-1 (c)∈R 49 Directional consistency in the spatial domain. The cosine similarity between two spatial feature maps ranges from -1 to 1. Values ​​closer to 1 indicate higher similarity between images in the feature space. Conversely, values ​​closer to -1 indicate lower similarity between images in the feature space. Figure 2 (e) and (f), we can see that the spatial feature map F s，t (c), F s，t-1 (c)∈R 49 The similarity between them is too high, such as Figure 2 As shown in (f), the similarity feature S t-1，t tends to 1 in all dimensions, which will reduce the effectiveness of subsequent regression. Figure 2 (e), differential feature D t-1,t It has more significant feature dimensions, which provides a better basis for subsequent regression.

[0163] 2.5) Using the fully connected regression module, regress and fuse the difference features obtained in step 2.3) and the similarity features obtained in step 2.4) to obtain the trained odometry estimation model (GPR-OdomNet);

[0164] The fully connected regression module integrates two-way features for motion prediction, including parallel regression and adaptive fusion operations. First, the above difference feature D t-1,t and similarity feature S t-1，t Independent regression, the formula is as follows:

[0165]

[0166] Among them, W d ,W s ∈R 512 is the learnable weight, b d , b s is the learnable bias value; is used to quantify the absolute intensity of feature changes, It is used to measure the consistency of feature direction; finally, the above absolute strength is adjusted by adaptive weight α and consistency Perform dynamic fusion to obtain the travel distance O of the vehicle or robot t-1，t , the formula is as follows:

[0167]

[0168] The adaptive weight α is used to adjust the influence of the difference and similarity features on the distance estimation. Finally, the training is completed to obtain the trained odometry estimation model.

[0169] 3) Using a ground-penetrating radar (GPR) mounted on a moving vehicle or robot, a series of GPR raw signals of underground features are collected. These signals are then processed according to the method in step 1) to obtain multiple B-scan images. Two adjacent B-scan images are then fed into the trained odometry estimation model to output the distance traveled by the vehicle or robot.

[0170] In order to verify the effect of the present invention, the inventors conducted the following experiments:

[0171] The odometry estimation model (GPR-OdomNet) employed in this paper utilizes both the differences and similarities of B-scan images for distance regression, significantly improving the accuracy of Euclidean distance estimation between adjacent frames. Compared to the commonly used learning-based GPR odometry model, it achieves a 9.1% reduction in root mean square error (RMSE), as shown in Table 1.

[0172] Table 1 RMSE (cm) of different relative distance estimation methods based on GPR

[0173]

[0174] The original GPR signals in the CMU-GPR dataset used in this paper come from basement / factory / parking lot scenes and are implemented on the NVIDIA RTX 3090 platform.

[0175] Since GPR odometry is often used in combination with other sensors such as IMUs and wheel encoders, it is necessary to compare localization performance in a sensor fusion setting. This paper compares the performance of GPR-OdomNet and uses a factor graph model fused with IMUs and wheel encoders as a baseline method. Then, as Figure 3 As shown in Figure 3, GPR-OdomNet is incorporated into the factor graph model. Within the adopted factor graph optimization framework, the IMU observations are processed using a pre-integration method to constrain the robot's position, orientation, and velocity, while the measurements from the wheel encoders and the estimates from GPR-OdomNet constrain the robot's velocity.

[0176] exist Figure 3 In the figure, square nodes represent factor nodes, and circular nodes represent variable nodes. The Prior (PR) node is a prior factor node. The IMU (IM) node primarily constrains the attitude and velocity between two consecutive moments. The Wheel Encoder (WE) and GPR-OdomNet (GO) nodes constrain the velocity at the next moment t. The Offset (BI) node constrains the offset between two consecutive moments.

[0177] The comparison results are as follows Figure 4 As shown in Figure 2, the ground truth data refers to the data collected by the total station. The method without GPR-OdomNet refers to the baseline method. The method with GPR-OdomNet refers to the method that combines IMU, wheel encoders, and GPR-OdomNet. Figure 4 (a) and (b) are two trajectories from the gates_g dataset: Figure 4 (a) is a trajectory that cycles twice, Figure 4 (b) is a trajectory that first goes back and forth vertically, then moves horizontally, and finally moves vertically again. Figure 4 (c) is the round-trip trajectory of the nsh_b dataset. Figure 4 (d) is a slightly curved back-and-forth trajectory from the nsh_h dataset. The results show that the methods with GPR-OdomNet significantly outperform those without GPR-OdomNet.

[0178] The results show that the GPR-OdomNet method significantly improves the odometry accuracy on all three datasets (gates_g, nsh_B, and nsh_h).

[0179] To further compare the odometry performance, Figure 5 Shows the Figure 4 The RMSE value of the absolute trajectory error (ATE) corresponding to the trajectory in is:

[0180]

[0181] Where T represents the total duration of the movement, x t and y t denote the horizontal and vertical coordinates obtained from the total station at time t, respectively, and and Then denote the horizontal and vertical coordinates estimated at time t, respectively. Figure 5 (a)-(d) correspond to Figure 4 Trajectories in (a)-(d). The results show that the model with the GPR module achieves significantly lower error than the model without it, and this difference becomes more pronounced over time. This metric provides a global assessment of localization and ranging results, making it more comprehensive. Clearly, the GPR-OdomNet approach outperforms the other methods. Furthermore, Table 2 shows that the GPR-OdomNet approach achieves the best accuracy.

[0182] Table 2 Comparison results of ATE RMSE(m) for each data set

[0183]

Claims

1. A method for estimating odometer based on ground penetrating radar, characterized by: The method for estimating the odometer of a ground penetrating radar comprises the following steps performed in sequence: 1) Using the publicly available CMU-GPR dataset, the GPR raw signals were subjected to background removal, DC offset removal, triangular bandpass filtering, SEC gain compensation, wavelet denoising, and Gaussian smoothing to obtain multi-frame B-scan images; 2) constructing an original odometry estimation model, and then inputting any two adjacent B-scan images obtained in step 1) into the original odometry estimation model for training to obtain a trained odometry estimation model; 3) Using a ground-penetrating radar (GPR) mounted on a moving vehicle or robot, a series of GPR raw signals of underground features are collected. These signals are then processed according to the method in step 1) to obtain multiple B-scan images. Two adjacent B-scan images are then fed into the trained odometry estimation model to output the distance traveled by the vehicle or robot.

2. The method for estimating odometer based on ground penetrating radar according to claim 1, characterized in that: In step 1), the method of performing background removal, DC offset removal, triangular bandpass filtering, SEC gain compensation, wavelet denoising and Gaussian smoothing on the GPR original signal to obtain a multi-frame B-scan image is: 1.1) Background removal Calculate the average value of all scan signals on the horizontal track and subtract it from the GPR original signal to obtain the background removal component. The formula is as follows: Where N is the number of horizontal channels, t is the time sampling point, x is the horizontal channel number, and s(t, x) is the xth row of the original signal; 1.2) DC offset removal Perform cubic polynomial curve fitting on each scan trace to obtain the baseline component: The baseline component is then subtracted from the GPR original signal s(t) to obtain the component after removing the DC offset, as shown in the following formula: s dewow (t)=s(t)-s base (t) Among them, c k are the polynomial coefficients; 1.3) Triangular bandpass filtering The filter coefficients are determined by the following equation: Among them, f min =200MHz, f max =850MHz is the passband boundary, f s is the sampling frequency, w tri is a triangular window function, n is an integer variable, and its value range is -K≤n≤K, where K=(N-1) / 2, and N is the filter length or order; Based on the filter coefficients and the GPR original signal s(t), the filtered components are calculated as follows: s filtered (t)=s(t)*h[n] 1.4)SEC gain compensation processing First, define the time-varying gain function: G(t)=t b ·e a·t Among them, a is the exponential gain coefficient; b is the power law gain coefficient; then the above time-varying gain function is used to calculate the value of the time threshold value t th The GPR original signal s(t) is gain saturated to obtain the time-varying gain function G after limiting clip (t): Based on the time-varying gain function G after limiting clip (t) and the GPR original signal s(t), and the enhanced component is calculated: s sec (t)=s(t)·G clip (t) 1.5) Wavelet denoising The GPR original signal s(t) is decomposed based on the Daubechies wavelet basis function, and the formula is as follows: c j，k =〈s,ψ j,k > The Bayesian adaptive threshold algorithm is used to determine the optimal threshold: in, is the noise variance estimate; Based on the Daubechies wavelet basis function and the optimal threshold, the wavelet denoising component is obtained: 1.6) Gaussian smoothing The Gaussian kernel function with a standard deviation of σ is used to perform spatial smoothing on the GPR original signal s(t). The formula is as follows: The time domain smoothing is achieved through convolution operation to obtain the smoothed component:

3. The method for estimating odometer based on ground penetrating radar according to claim 1, characterized in that: In step 2), the original odometer estimation model is constructed, and then any two adjacent frames of B-scan images obtained in step 1) are input into the original odometer estimation model for training. The method for obtaining the trained odometer estimation model is: 2.1) Constructing a raw odometry estimation model; the raw odometry estimation model consists of a cascade of a feature extraction module, a parallel difference detection module, a similarity detection module, and a fully connected regression module; wherein the feature extraction module consists of a cascade of a backbone network and an adaptive average pooling layer; the backbone network uses the first six layers of the ResNet-18 network from the initial convolutional layer to the fourth residual block; the difference detection module consists of a cascade of a tensor subtraction layer, a 3×3 convolutional layer, a batch normalization layer, an activation function, and a dropout unit; the similarity detection module consists of a cascade of a flattening layer and a cosine similarity calculation unit; and the fully connected regression module consists of two parallel fully connected layers; 2.2) Inputting any two adjacent B-scan images obtained in step 1) into the feature extraction module in the above-mentioned original odometry estimation model to extract features and obtain global features; First, the backbone network in the feature extraction module is used to extract the B-scan images B at time t-1 and time t through stacking operations of convolution, batch normalization, and ReLU activation. t-1 and B t The spatial features in the 1-channel 224×224 grayscale B-scan image are converted into a 512-channel 7×7 spatial feature map. The formula is as follows: F s =ReLU(BN(Conv 3×3 (I)))+Identity(I) Where I∈R B×1×224×224 is the batch of B-scan images, ReLU is the ReLU activation function, Identity(I) is the residual connection, Conv is the convolution operation, and BN is the batch normalization operation; Then, the adaptive average pooling layer is used to globally compress the above spatial feature map. The formula is as follows: Among them, H is the height of the spatial feature map, W is the width of the spatial feature map, and c is the channel index; Get the global feature F g ∈R B×512×1×1 ; 2.3) Using a difference detection module to perform difference detection on the global features obtained in step 2.2) above to obtain difference features; The tensor subtraction layer in the difference detection module is used to calculate two global features F g，t-1 and F g，t The absolute difference between them is used to obtain the difference vector representing the feature difference. The formula is as follows: ΔF=|F g，t -F g，t-1 | Based on the above difference vector ΔF, a 3×3 convolutional layer is used to extract the local difference features of the two global features. The formula is as follows: Where W∈R 512×512×3×3 is a learnable convolution kernel, represents the convolution operation, b c is the bias term; then based on the above local difference features, the batch normalization layer is used to stabilize the training process and obtain the normalized features. The formula is as follows: Among them, μ B ,σ B is the batch statistic, ε is a very small constant to prevent division by zero, and β is a learnable bias parameter; Finally, based on the above normalized features, we use the activation function and the discard unit to introduce nonlinear activation and regularization to obtain the difference features. The formula is as follows: D act =max(0.2×D norm ,D morm) D t-1,t =Dropout 0.2 (D act ) Among them, the 20% dropout rate Dropout improves the generalization of the model through random inactivation; 2.4) Using a similarity detection module to perform similarity detection on the spatial feature map obtained in step 2.2) above to obtain similarity features; The spatial feature map F obtained in step 2.2) is transformed into s Flattened, the formula is as follows: f:R B×512×7×7 →R B×512×49 Then, the cosine similarity calculation unit is used to calculate the spatial feature map F channel by channel. s The cosine similarity of , continues to obtain the similarity feature, the formula is as follows: 2.5) Using a fully connected regression module, regress and fuse the difference features obtained in step 2.3) and the similarity features obtained in step 2.4) to obtain a trained odometry estimation model; First, the above differential features D t-1，t and similarity feature S t-1，t Independent regression, the formula is as follows: Among them, W d , W s ∈R 512 is the learnable weight, b d , b s is the learnable bias value; is used to quantify the absolute intensity of feature changes, It is used to measure the consistency of feature direction; finally, the above absolute strength is adjusted by adaptive weight α and consistency Perform dynamic fusion to obtain the travel distance O of the vehicle or robot t-1，t , the formula is as follows: Finally, the training is completed to obtain the trained odometry estimation model.

Images