Research on Locating Tunnel-Lining Defects Using Fast Synthetic Aperture Focusing Imaging Based on GPR

. This paper proposes a method based on image processing algorithms for ground penetrating radar (GPR) to locate hidden defects in tunnel linings. Firstly, the fast synthetic aperture focusing imaging (Fast-SAFI) algorithm is used to accurately identify the morphology of tunnel-lining defects. Secondly, an iterative algorithm is used to determine the connected regions on the binary image, exclude background noise interference, and locate the centroid and vertices of the correct target connected regions to achieve the positioning of the depth of tunnel-lining defects. To verify the feasibility of the proposed positioning algorithm, a verification experiment was conducted on the experimental wall of the China Academy of Railway Sciences. The experimental results show that the proposed positioning algorithm is reliable and rapid for identifying and locating the morphology and depth of tunnel-lining defects


Introduction
Tunnel lining is an important component of tunnel engineering.However, it is often subject to geological conditions and environmental influences, resulting in various hidden defects such as cracks, voids, and disconnections.The hidden defects in the tunnel lining not only cause serious damage and instability to the structure of tunnel linings, but may also cause serious accidents.Therefore, discovering and solving these potential problems is crucial for the safe operation of tunnels [1], [2].Ground penetrating radar (GPR) is an efficient, anti-interference, and high-penetration tool for non-destructive testing of underground structures, providing more guarantees for the safe and scientific construction of underground engineering [3].By transmitting electromag-netic waves and receiving reflection signals, GPR can form B-scan profiles to interpret and analyze tunnel-lining defects, and infer their location, morphology, and type [4], [5].However, tunnel-lining defects are complex in shape and diverse in type, and relying solely on manual interpretation is difficult to accurately identify and locate defects, leading to problems such as missed diagnosis and misdiagnosis [6].In order to improve tunnel safety and reduce inspection costs, automated detection has become an important means of infrastructure inspection and is gradually becoming a trend for future development [7], [8].
Synthetic Aperture Radar (SAR) is a radar-based imaging technology that can achieve high-resolution imaging of ground and ocean targets.The SAR uses the synthetic aperture generated by its own motion to form a large aperture equivalent to the actual aperture, thus improving the imaging resolution.Through the SAR technology, high-resolution and high-precision information of ground objects can be obtained, and it is widely used in military, civilian, aerospace, and other fields.Similar to the SAR radar system, the GPR is also based on electromagnetic principles and radar equations.Usually, the GPR systems operate in a single-point transmitreceive mode, so the synthetic aperture focusing technique (SAFT) can also be applied to GPR [9].Because the wave equations for elastic waves and electromagnetic waves are similar, the GPR synthetic aperture imaging algorithm is usually evolved from seismic imaging algorithms [10].The purpose and requirements of these two imaging methods are also very similar.Currently, the main types of migration algorithms for the GPR include frequency-wavenumber migration, phase-shift migration, Kirchhoff migration, and backpropagation (BP) algorithm.Frequency-wavenumber migration uses Fast Fourier Transform (FFT) and is suitable for fast focusing of single-layer homogeneous media flaws [11].Phase-shift migration algorithm has high complexity and requires very accurate parameters for wave velocity and layer thickness, which is not suitable for processing actual ground penetrating radar data [12].Kirchhoff migration generates a large amount of noise during the implementation of GPR B-scan focusing, which can easily mask the reflected echoes of underground flaws [13].Back-projection-based focusing method is a time-domain migration algorithm that uses geometric equations for implementation.Compared with other imaging algorithms, its imaging model is relatively simple and has strong applicability [14].
The GPR migration algorithms can achieve good imaging results in practical applications.For example, Bradford [15] proposed a reverse-time migration (RTM) algorithm that can handle rugged terrain and varying surface heights to generate high-fidelity underground structural images, which has been applied in dune exploration.Giannakis, et al. [16] developed an improved Kirchhoff migration algorithm to evaluate the internal structure of tree trunks and detect the degree of decay for locating internal tree diseases.Alani [17] utilized the inverse-time migration imaging technique to accurately locate diseased areas in the mapped tree trunk, and the effectiveness of the imaging technique was demonstrated through simulation and field experiments.Tanyer [18] proposed a delay-weighted sum algorithm to study the effect of non-uniform media on the inversion BP imaging, which effectively improves the imaging quality.
In recent years, the surge in deep learning has introduced new approaches to image processing, and Convolutional Neural Networks (CNN) have become one of the methods for processing GPR images and performing defect recognition.For instance, ResNet-50 [19] is adopted to detect the tunnel lining structural condition.Rosso, et al. [20] compare three deep learning models: ResNet-50, EfficientNet-B0, and ViT.And the detection accuracy was compared with and without preprocessing.However, an effective localization strategy is yet to be proposed.An algorithm combining Self-Attention Dense Contrastive Learning (SA-DenseCL) and Mask Region-Convolutional Neural Network (Mask R-CNN) [21] is also proposed, capable of defect localization after training.
The existing methods can accurately identify and locate the morphology and depth of reinforcement but lack a precise identification and localization method for complex and irregular disease morphology and depth.To solve the problem of difficult identification and positioning of the tunnel-lining disease morphology and burial depth, this paper proposes a fast synthetic aperture focusing algorithm to determine the disease morphology, combined with image binarization and target vertex method to determine the disease burial depth.The effectiveness of the proposed method is validated by focusing imaging indicators on simulated and field data.This paper will introduce the principles of Synthetic Aperture Focusing Imaging, as well as an enhanced Fast Synthetic Aperture Focusing Imaging method that is more effective and suitable for synthetic aperture data in Sec. 2. Subsequently, the localization method will be discussed in Sec. 3, encompassing image binarization applied to the images produced by the Fast-SAFI technique, along with the target vertex method used for locating the burial depth of the defects.Section 4 discusses the reliability and complexity of the algorithm.

Synthetic Aperture Focusing Imaging
In the scanning process of the GPR, the detected electromagnetic waves will generate different travel times during propagation, causing any scatterer in the B-scan image to be displayed as a hyperbola.Due to the long tail of the hyperbola, the resolution along the synthetic aperture direction mainly depicts low-resolution features.Therefore, it is challenging to locate defects in tunnel linings using unfocused B-scan images, and it is necessary to transform the scattered hyperbolas into a focused image that displays the actual position, size, and corresponding electromagnetic wave reflectivity of underground objects.Therefore, the GPR imagefocusing techniques are commonly used to fold the detected hyperbolas into a concentrated spot at the vertex, highlighting the underground targets.
The fundamental principle of the synthetic aperture focusing technique in the GPR is to transform the imaging of underground targets into the focusing of the targets.Hyperbolas can be focused into smaller points, achieving precise positioning of underground targets.This can improve the azimuthal resolution, weaken interference, and provide accurate data.By using the focusing technique, A synthetic aperture focusing imaging scene model for the GPR data has been constructed, as shown in Fig. 1.In this scene model, the detection direction is along the x-axis, and the depth direction is along the y-axis.The GPR transmitter is at a distance of ℎ from the lining surface and samples uniformly along the detection direction to obtain GPR data   (  , ).
Assuming the GPR transmitter is located at point  (  , −ℎ) and the electromagnetic waves are reflected at point  ( P ,  P ) and refracted at point  ( R , 0) before being received back at point , the distance traveled by the electromagnetic waves from  to  can be calculated as follows: When the focal point  is located at point , the delay of data   (  , ) can be expressed as follows: where  represents the speed of light, and  represents the velocity of electromagnetic waves within the lining.According to Snell's law, the incident angle  1 and the refracted angle  2 can be expressed as follows: sin According to (4) and considering coordinate transformations, a fourth-order polynomial equation with respect to the horizontal coordinate  P of the focusing point  can be derived as follows: The abscissa  P of the focusing point , which can be obtained based on (5), is substituted into (2) to obtain the time delay  P .
After collecting  data points at each focal point from the aperture, the amplitude of the reconstructed image points can be obtained using the concept of "delay-and-sum".This algorithm adds the signals at the focal points, effectively combining the signals from multiple scattered points and improving the imaging quality.

Fast Synthetic Aperture Focusing Imaging
Due to the larger number of aperture channels and different baseline channel counts of the focusing points involved in the computation below the image, the data cannot be processed further.To improve the applicability of synthetic aperture data, a fast synthetic aperture focusing Imaging (Fast-SAFI) algorithm based on delay-and-sum averaging is proposed to address the issue of varying baseline channel counts during data processing.This greatly improves the efficiency of GPR synthetic aperture focusing imaging.The specific improvement method is as follows: 1. Simplify the calculation of reflection points using the Johansson approximation formula [22].This simplifies the reflection point formula into a linear judgment formula, where only half of the reflection points and delays at each focal point need to be calculated, and the other half can be obtained through horizontal symmetry to reduce computational costs.
2. Simplify the calculation of the delay matrix.The GPR data is represented by a matrix [ × ], with each channel containing sampling points as focusing points and a channel spacing of Δ.Assuming that a focusing point on a certain channel is  ( P ,  P ), the number of channels in the synthetic aperture is , the reflection point coordinate matrix is Q P (  , ℎ), and the delay matrix is n ′ .In addition, focusing points at the same depth have the same number of reflection points Q P (  + Δ, ℎ), and the delay matrix n ′ is also the same.Let the delay matrix of one channel be the same as that of the other  − 1 channels.Therefore, only three parameters, namely the distance between the transmitter and the surface h, the number of sampling points, and the dielectric constant of the lining, are required to calculate the delay matrix for each channel.
3. Encapsulate the delay matrix groups for all channels into a "cell" and store it as a whole on the computer, avoiding the need to calculate the delay matrix for each focal point.This greatly reduces the computational workload and significantly improves the calculation speed, making it suitable for GPR data processing.
The input for the fast synthetic aperture focusing technique (Fast-SAFI) algorithm is the GPR B-scan matrix s = [ × ], with an inter-scan distance of Δ and a time sampling interval of Δ.The transmitter height above the interface is ℎ.The sampling interval in the subsurface layer is Δ, and the relative dielectric constant of the subsurface layer is .The flowchart of the Fast-SAFI algorithm is shown in Fig. 2.
Step 1: Calculate the delay matrix for all columns of GPR data, looping through  columns.
Step 2: Calculate the number of synthetic aperture channels _number + 1 and the synthetic aperture length .
Step 3: Determine if the number of channels at the point is 1.
If the number of channels is 1, the delay matrix at the point is set to  p .Otherwise, set _number/2 as the number of times to move the data channel above the focusing point to the left and right, looping through _number/2.Step 4: Using _number/2 values, calculate the reflection point and assign values to the delay matrix from the center to both ends, obtaining the delay matrix  p .
Step 5: Package the delay matrix as an "A cell" and store it in the computer.
Step 6: Read in the other  − 1 data channels of the GPR B-scan matrix s = [ × ], loop through rows, and use "A cell" for looping.
Step 7: Find the corresponding delay matrix for the focusing point and find the values of each synthetic aperture channel .
Step 8: Sum the delayed matrix elements of  and take the average to obtain the imaging result of the focusing point.
Step 9: Inversion imaging of GPR defects.

Binarization
After applying the Fast-SAFI technique, the pixel energy distributed on the hyperbolic curve can be summed up at the real location of the subsurface anomaly.However, when the image is shifted, background information still causes significant interference to the target, so it is necessary to further enhance the difference between the target and the background.In this regard, image binarization methods can be used to segment the target and the background for more accurate target positioning.This step can effectively reduce the impact of noise, improve the accuracy of subsurface anomaly location, and facilitate further research and analysis.
Image binarization is an image processing technique that simplifies and separates the image, and there are various threshold selection methods, including bimodal method, maximum inter-class variance method, maximum entropy thresholding, and iterative method.Each method has its unique advantages and applicable scope, and the most suitable method should be selected based on actual situations for image binarization.
The bimodal method [23] begins by converting the input image to grayscale.It then constructs a histogram of pixel intensities.The algorithm identifies two peaks in the histogram, representing the foreground and background intensities, and selects a threshold value that separates these two modes.Maximum inter-class variance method [24] computes the histogram of the input grayscale image.It iteratively evaluates the variance between two classes (foreground and background) for different threshold values and selects the threshold that maximizes this inter-class variance.Maximum entropy thresholding [25] calculates the entropy of the binary image for different threshold values.It seeks the threshold that maximizes the entropy, indicating the highest information content.This balance between object and background regions is achieved by maximizing the uncertainty in the image.
The basic idea of the iterative method is to first select an initial threshold value, and then use this threshold to separate the image into two parts: the target and the background.The average grayscale values of these two parts are then calculated.The average of these two values is used as the new threshold, and the process is repeated until the new threshold is equal to or only slightly different from the previous iteration's threshold.In other words, the convergence condition of the iteration is met.The advantage of the iterative method is that it can adaptively select thresholds for images with different grayscale distributions.Moreover, since this method only relies on the statistical information of the image itself and does not require human intervention, it has good robustness.The flowchart of the iterative threshold algorithm is shown in Fig. 3, and the specific steps of the process are as follows: Step 1: Initialization.Set the threshold value  0 for all pixels in the focused GPR B-scan to an initial value,  0 = 0.6 × ( max +  min ), where  max represents the maximum gray value in the image, while  min represents the minimum gray value in the image.
Step 2: Segmentation.Compare the grayscale value of each pixel in the focused GPR B-scan with the initial threshold value  0 .If the grayscale value is less than  0 , mark it as background; otherwise, mark it as foreground.
Step 3: Update the threshold value.Calculate the average grayscale values of the foreground and background pixels, and use the average of these values as the new threshold  h .
Step 4: Convergence judgment.Compare the difference between the new threshold value  h and the previous threshold value  0 .If the difference is less than the predetermined threshold, the algorithm converges and outputs the final segmentation result.Otherwise, return to Step 3 and continue the iteration.
Step 5: Input the GPR B-scan after iterative thresholding.The final result is a black and white binary image, which is convenient for subsequent processing and analysis.
In order to discuss the convergence of the iterative method, the grayscale values of all  pixels can be rearranged into a sorted list in ascending order { 1 ,  2 , . . .,   }.By comparing each value with  0 , the list can be divided into two parts.The average values of these two parts   1 ,   1 can then be computed: (7) where  1 represents the number of pixel points with grayscale values less than the threshold  0 .Then, the new threshold  h , referred to as  1 , can be obtained: Obtaining the threshold   after  iterations is straightforward: where   represents the number of pixel points with grayscale values lower than the threshold  −1 .
The convergence criteria |  −  −1 | can be expanded as follows: The experiments indicate that when the initial value is set as  0 = 0.6 × ( max +  min ), for the B-Scan images used in this experiment,  1 >  2 > • • • >   , and it can be observed that which means that the sequence |  −  −1 | converges.Therefore, the iterative method used to binarize images in this experiment is proven to be feasible.
After comparing the four different threshold selection methods mentioned above, the results are shown in Fig. 4. It was found that the iterative method had the best performance in binarizing the four offset images.The binarized image processed by the iterative method presented white areas similar to the focused cluster of the original image, while the background outside the white areas was very clean with no obvious small color blocks.On the contrary, the other three threshold methods yielded poor results, with the target and background not completely separated and many large color blocks appearing, which may lead to a misjudgment in choosing the target.Specifically, the results obtained using the Maximum Entropy thresholding method exhibit significant information loss regarding the target vertices, leading to a considerable localization error.Therefore, the iterative method was chosen as the thresholding method for image binarization.

Depth Estimation
After binary thresholding the image, the target region is presented as white, while the background is black.To accurately locate the target, it is necessary to obtain the coordinate information of the white region.For this purpose, the connected regions need to be determined, the centroid position needs to be calculated, and the buried depth of the underground disease can be estimated.
In order to distinguish between the black and white regions in the binary image, they need to be labeled.However, when using the iterative method to process the image, there may be small and disconnected white regions that affect the selection of the target position.To solve this problem, all independent white connected regions were labeled, and a different index value was assigned to each closed white connected region for differentiation.After determining the target connected region, the centroid of this region needs to be calculated.By utilizing the coordinates of all pixels, the centroid position and height information of the bright cluster can be obtained, which helps to better understand the shape and features of the target region.
The horizontal and vertical coordinates of the centroid are shown below: where  represents the total number of pixels in the white connected region,   represents the horizontal coordinate of pixel , and   represents the vertical coordinate of pixel .
By calculating the minimum vertical coordinate  min and the maximum vertical coordinate  max of the white connected region, the height  of the white connected region is obtained: The vertex of the white connected region is located at a height /2 above the centroid.The expression for the coordinates of the vertex position is shown below: After obtaining the coordinates of the vertex of the white connected region, the vertex is fitted to the extracted hyperbolic curve target to observe the accuracy of the vertex positioning.Figure 5 shows the results of each processing step of the proposed localization algorithm.The results show that the proposed algorithm can accurately locate the red localization point on the vertex of the hyperbolic curve, effectively solving the problem of vertex localization of tunnel-lining diseases.

Performance Evaluation
To evaluate the performance of algorithms, parameters such as Integrated Sidelobe Ratio (ISLR), Signal to Clutter Ratio (SCR), and computational speed are integrated.
The ISLR is defined as the ratio of the total energy within the sidelobes to the peak energy of the main lobe: The SCR is also a key parameter for evaluating the quality of focusing ability.It can be defined as the ratio of the received target signal power to the received clutter power, as shown below: where  target is the received target signal power, and  total is the total power of the image.

Pattern Recognition Results
In order to verify the effectiveness of the Fast-SAFI imaging algorithm, tunnel-lining models were established and the BP algorithm was used as a comparative algorithm to analyze the focusing imaging effects of three types of diseases: cracks, voids, and delamination.Tunnel-lining disease models and their corresponding forward images are shown in Fig. 6.For cracks, voids, and delamination diseases, both the BP algorithm and the Fast-SAFI algorithm can improve the visibility and recognition of the diseases.The results of the focusing algorithm are shown in Fig. 7.
As shown in Fig. 7(a), the BP algorithm can refine and clarify the edge and texture features of the defects, making them more visible and recognizable.In addition, the BP algorithm can optimize and adjust the background around the defects to highlight their contours and features, thereby improving the identification of the defects.However, the BP imaging algorithm may face some difficulties in handling complex voids, and delamination defects.This is because the shape and size of the voids, and delamination defects are irregular and usually accompanied by local color and brightness changes.These problems may cause some distortion and misjudgment in the BP imaging algorithm when dealing with voids, and delamination defects.
As shown in Fig. 7(b), the Fast-SAFI algorithm has certain advantages over the BP algorithm in processing cavity and void defects.Void defects are focused into an approximate circular shape, which can better estimate their burial depth.Delamination defects are difficult to identify their approximate shape due to significant fluctuations in the surrounding dielectric constants caused by their proximity to the surrounding rocks.However, after the Fast-SAFI algorithm, the reflection signal on the surface of Delamination defects is clear, which can determine the approximate burial depth and depth range.Table 1 summarizes the processing time, Signal-to-Clutter Ratio (SCR), and Integrated Side Lobe Ratio (ISLR) performance of the two proposed algorithms in simulated tunnel-lining models.Compared with the BP algorithm, the Fast-SAFI algorithm has faster computation time, strong clutter suppression ability, and excellent anti-interference ability.
Table 2 illustrates the comparison between the Fast-SAFI algorithm and existing B-Scan imaging methods, such as HSA, PSA, -kA, and BPA [26].It can be observed that Fast-SAFI outperforms the existing methods in terms of the SCR metric.

Localization Results
Firstly, the Fast-SAFI algorithm is used to determine the true location and shape of the underground defects.Next, image binarization technique is applied to determine the target connected area.By calculating the coordinates of all pixels, the centroid position, vertex position, and height information of the defects can be obtained.Ultimately, the automatic localization of the shape and burial depth of tunnel-lining defects is achieved.
The complexities of the four mentioned image binarization methods above, namely, bimodal method,maximum inter-class variance method, maximum entropy thresholding, and iterative method are all  ().Now represent them individually as  (  ()),  (()),  (ℎ()), and  (()) to compare their complexities, where  represents the side length of the input image where  represents the number of iterations.
Therefore, it can be observed that, within a certain range of input image dimensions, the complexity of the bimodal method and the iterative method with fewer iterations is more competitive than the complexity of the maximum inter-class variance method and the maximum entropy thresholding.
When processing images with an input size of 465×465 using MATLAB 2020b on a computer with a 12th Gen Intel®Core™i5-12490F CPU, although the iterative method is more time-consuming, it is still within the same order of magnitude as the other three methods.Table 3 compiles the average processing time for individual images when using the four mentioned binarization methods on the same set of GPR B-Scan images.In this project, the use of the iterative method as the image binarization technique, although more time-consuming compared to the other three methods, has minimal impact on the system's overall performance due to its small computational scale.However, it offers superior target detection capabilities and background noise suppression, ultimately enhancing the accuracy of defect localization.Therefore, it appears to be a feasible choice.
In order to validate the feasibility of the proposed algorithm for identifying the morphology and locating the depth of tunnel-lining defects, field experiments were conducted using the railway institute's experimental wall as the test object.In the GPR B-scan image, a total of four target prediction boxes were generated.Among them, the hyperbolic target features in boxes 1, 2 and 3 were very clear and could be easily identified with the naked eye.However, the hyperbolic features in the selected area of box 4 was not very obvious, making it difficult to determine whether there was a defect target through manual judgment.Nonetheless, under the assumption that box 4 was not a false detection, the proposed localization algorithm was used to locate the defect targets in the regions selected by these four boxes, as shown in Fig. 8.In the Fast-SAFI images of boxes 1, 2, and 3, irregular diseases are clearly restored into circular targets, which can effectively distinguish the target from the background.The results of the localization algorithm for the four target regions are shown in Fig. 9.
However, in box 4, there are no obvious hyperbolic features and it is difficult to distinguish between the target and the background.This may be because there are interference anomalies such as broken stones in the area where boxes 2 and 3 are located, and there are diseases in the area.After the Fast-SAFI algorithm is processed, multiple clutter energy points are formed, which further increases the possibility of misidentification in box 4. The focusing result can also serve as a reference for judging the existence of underground diseases.After processing with the localization algorithm, the exact positions of each target within their respective bounding boxes have been determined.Figure 10 shows the localization results of the regression of the original image, with the target vertices of each bounding box marked as solid red dots.Overall, including one suspected false positive, all targets had a depth of greater than or equal to 0.5 meter.
In theory, this system can be deployed as a real-time system on GPR vehicles, the minimum driving speed compatible with the system can be expressed by the following equation: where  is the detection length per image,  is the maximum average time required for imaging and detecting defects,  is the number of defects, and  is the maximum average time for detecting defect locations.Since in practical scenarios,  is very small, the term  is significantly smaller than .Therefore, the above equation can be approximated as The system proposed in this paper theoretically can be integrated as a real-time system and operate on GPR vehicles traveling at speeds of approximately 20 km/h.

Conclusions
In this paper, we have proposed a method for identifying the morphology and depth of tunnel-lining defects for GPR data.The method uses the Fast-SAFI algorithm to achieve high-resolution focused imaging, and iteratively estimates the vertex position of the defect using target localization estimation method, thereby obtaining the true shape and depth of the defect.Although using the iterative method for image binarization incurs higher computational time, from the perspective of the entire system, this additional time can be almost negligible.Moreover, the benefits it brings are clearly evident.The processing results from the field measured data show that our method is more stable and convenient, which can meet engineering requirements.Through the application of this method, the problem of inaccurate positioning of tunnel-lining defects in shape and depth can be effectively solved, providing a feasible solution for engineering practice.
Our work still has some limitations, and in the future, we will continue to improve the accuracy of localization and the lightweight nature of the algorithm.We will also explore the integration of relevant deep learning models to further optimize the system into a real-time solution, better meeting the requirements of practical detection.

Fig. 1 .
Fig. 1.A synthetic aperture focusing imaging scene model for the GPR data.

Fig. 2 .
Fig. 2. The flowchart of the Fast-SAFI algorithm.As a reference, the input GPR B-scan matrix s[, ] corresponds to the image shown in Fig. 6(b), and the output GPR B-scan data I[, ] corresponds to the image shown in Fig. 7(b).

Fig. 9 .
Fig. 9.The results of the localization algorithm for the four target regions.

Fig. 10 .
Fig. 10.The localization results of the regression of the original image.
Performance of the focusing algorithm on different models.Comparison of SCR for different methods.