Automated Defect Detection and Visualization for the Robotic Airport Runway Inspection

Detection of both surface and subsurface defects is a vital task for maintaining the structural health and reliability of airport runways. We report the automated data collection and analysis for airport runways based on our novel robotic system, which employs a camera and a GPR (Ground Penetrating Radar) to inspect the surface and subsurface conditions, respectively. To perform the automated data analysis, we propose a novel crack detection algorithm based on the images, and a subsurface defect detection method with GPR data. Additionally, to create a composite global view of a large airport runway span, a camera/GPR data sequence from the robot is aligned accurately to create a continuous mosaic for visualization. We combine these algorithms into a software to perform automated on-site analysis. We have put our robot and software into engineering practice over 20 airports in China, achieving the performance of 70% and 67% F1-measure for crack detection and subsurface defect detection, respectively. More importantly, the results of our algorithms can satisfy the requirement of applications.


I. INTRODUCTION
Structural defect inspection is essential task to monitor the health and reliability of the airport runways. Early detection of defects, such as surface cracks and subsurface voids, for airport runway is an important maintenance task. By the end of 2019, there are 238 civil transport airports in China, with nearly 300 runways. Structural defects exist in many airport runways, leading to costly maintenance even safety risks if serious defects cannot be detected in their early stages.
Traditional manual defect detection methods are very labor-intensive, time-consuming, with low-accuracy, and error prone. Even slight defects may be early warning signs of significant failure for airport runways, and need to be detected accurately and timely. Thus, it is necessary to develop an automated defect detection method for airport runway inspection.
We develop an airport runway inspection robot (ARIR), with a camera and a GPR (Ground Penetrating Radar) fixed to perform the condition sensing task. This paper focuses on addressing the problem of analyzing the collected data The associate editor coordinating the review of this manuscript and approving it for publication was Yizhai Zhang. accurately and automatically. Specifically, we propose a robust image-based crack detection algorithm, and a deep learning based subsurface defect detection method from GPR data. In addition, to create a composite global view of a large airport runway span, a camera/GPR data stitching algorithm is presented to create a continuous mosaic for visualization.
The rest of paper is organized as follows: the related work is summarized in Section II before we introduce our sensing suite configuration and data collection scheme in Section III. The proposed algorithms are illustrated in details in Section IV. The performance of our algorithms is validated in experiments in Section V, and the paper is concluded in Section VI.

II. RELATED WORK
In recent years, research in structure inspection using robotic systems has attracted much attention and resulted in several prototypes, such as the bridge deck inspection robot [1], [2], and the robot for tunnel structure health inspection [3]. However, the scenario of airport runway inspection is quite different and more challenging. We develop a robotic system for airport runway inspection. In order to implement the automated data analysis for airport runway inspection, there are three challenging problems, including large scale image stitching, crack detection from images, and defect detection from GPR data, that need to be addressed.
General image stitching techniques have been well studied [4]- [6] and even put into commercial use. The homography model is adopted in these methods and point features in images are used to estimate the homography model [7], [8]. However, the image stitching problem in airport runway inspection is significantly different in two points: On the one hand, the surveyed regions is very large, involving huge amount of images to be stitched. On the other hand, images of airport runway are usually lack of features, leading to occasional but inevitable failures of feature-based homography estimation.
Existing image-based crack detection methods can be mainly classified into two categories: the standard image processing methods, and machine learning based techniques. Standard image processing methods, including intensity thresholding [9], edge detection [10], and morphological filtering [11], [12], have been widely studied for automated crack detection. However, the performance of these methods is usually dependent upon the parameter choices [13] which are very difficult to accomplish for field images with significant visual clutters, leading to unreliable detection results in airport runway inspection applications. Machine learning based crack detection methods build on techniques such as support vector machines [14], random forest [15], random structured forest [16], and neural networks [17]. The machine learning based methods obtain more robust performance compared with image processing techniques. However, a supervised training stage is needed which requires large amount of labeled data, a difficult requirement for applications in airport runway since the sample images with cracks are rare. Our proposed crack detection algorithm combines the advantages of image processing methods in few samples and machine learning based techniques in robustness, achieving a satisfying performance on airport runway data.
Recognizing the subsurface objects automatically from GPR data is nontrivial, because a GPR cannot provide 3D positions but a reflection image full of significant signal clutter. Thus, object detection from GPR in an automated manner is still a challenging problem. Standard signal processing methods, such as template matching [18], S-transform [19], and wavelet transform [20], have been used in GPR data analysis. However, these methods are sensitive to noise, so cannot be employed directly in subsurface defect detection. Machine learning based methods, especially the CNN-based deep learning techniques [21], [22] obtain much attention in GPR data analysis [23], [24]. Although these CNN-based techniques have achieved a certain of positive results, they still cannot satisfy the requirement of field applications. One of the main reasons is that only 2D B-scan images are employed, which ignores the natural 3D property of subsurface defects.
Difference with all of the above mentioned works, our paper focuses on the airport runway data analysis which is collected by our ARIR system integrated with camera and GPR. The developed data analysis algorithms allows the robot to detect both surface and subsurface defects robustly, and build the large scale surveyed airport runway image to perform high-efficient assessments.

III. SENSOR CONFIGURATION AND DATA COLLECTION
This section firstly presents the sensor configuration of our ARIR system. The sensor setup of the ARIR system is shown in Fig. 1, which consists of one Raptor GPR with 900MHz antennas, and one DALS Nano M1920 camera with a resolution of 1920 × 1200.
To conduct defect inspection, the robot navigates within a predefined surveyed region on the airport runway to collect images and GPR data, as shown in Fig. 2. The robot first moves from the start to point A, then follows a linear path along each scan. The sensing suite can cover 1 meter wide in each scan. Once the robot finish the current scan, it moves to the next scan until the entire surveyed region is completely covered. At the end of each scan, the robot first moves to the turn point along an omni path, and then moves to the other scan to ensure there is no any region missed out. When scanning, the robot simultaneously transfers the image and GPR data to the nearby data analysis center in a van using 4G/5G connection. Then the collected data will be automatically analyzed off-line, which is presented in the following section.

A. LARGE SCALE IMAGE STITCHING FOR VISUALIZATION
There are two most challenging problem in our large scale image stitching. One is the large amount of images, easily causing the error accumulation; The other one is the lack of matched features between some adjacent images. To handle these problems, a two-stage image stitching algorithm is proposed. In this first stage, we align all images according to their residing GPS positions. Then in the second stage, we perform the adjustment under a global optimization using the available matched point features.

1) STAGE 1: POSITION BASED INITIAL STITCHING
Our robot has synchronized the RTK GPS and camera. Thus we can record global position for each image. Generally, the localization error of GPS while working in the open airport runway is about 5cm. Since we fix the camera downward to the ground, the image plane is approximately parallel to the horizontal airport runway plane. We build a world coordinate system {W } with X − Y plane to be horizontal. Then we align each image in the X − Y plane of {W } to form the initial stitching. Denote the resolution for each image as f x × f y pixels, and each image covers a region of w x × w y meters in the X − Y plane of {W }. Due to the GPS localization error, the initial stitching result need to be optimized in the following stage.

2) STAGE 2: FEATURE BASED REFINEMENT
In this stage, we first extract the matched point features between adjacent overlapped images. Then we define a cost function to perform global refinement.
Denote I a , a = 1, . . . , N I as the a-th image, with N I being the total number of images. We define I a and I b to be neighbors if their overlapping ratio is bigger than 0.2 in the original stitching result. Then we extract all matched points between I a and I b using SURF [25] method. The homography model is adopted to filter out the mismatched feature points under RANSAC (Random Sample Consensus) [26] framework. We define x a i and x b i to be homogeneous coordinates of two matched SURF points between two neighbored images I a and I b . A homography aims to map x b i to x a i following the relation where H a,b is the homography model between I a and I b . Thus the geometric distance between x a i to x b i can be computed as With the distance defined in (2), we find the outliers as the mismatched points iteratively under RANSAC framework. With the mismatched points removed, we keep the correctly matched SURF points for the following refinement.
In theory, the matched points from different images should be at the same position in the stitched image. Due to the errors of feature points, for each image I a , we combine all matched points with its neighbored images to compute the displacement for refinement. The displacement for I a is defined as where M a denotes the total number of I a 's neighbors, g(I a , j) is a function returning the image index number of the j-th neighbor of I a , S denotes the total number of matched feature points between I a and I g(I a ,j) , the symbol ↔ indicates one point corresponds to the other one. Thus, the global refinement for image stitching can be performed by solving the following optimization problem, arg min We summarize our larege scale image stitching algorithm in Algorithm 1.

B. CRACK DETECTION FROM 2D IMAGES
In a gray scale 2D pavement image, cracks can be visually distinguished from the rest of the image, because of their lower intensity values and their shape of continuous long curvilinear shapes. Based on these two characteristics, an automatic crack detection algorithm is proposed.

1) INTENSITY BASED CANDIDATE CRACK PIXEL EXTRACTION
First, the intensity property is used to extract the candidate crack pixels from the background. The image is first smoothed by a mean smoothing filter. We perform variations on the standard mean smoothing filter by means of threshold averaging, wherein smoothing is applied subject to the condition that the center pixel value is changed only if the difference between its original value and the average value is greater than a preset threshold. This has the effect that noise is smoothed with a less dramatic loss in crack like details. Find its neighbored images I g(I a ,j) ; 5 foreach I g(I a ,j) do 6 Extract corresponding points between I a and I g(I a ,j) using SURF; 7 Remove all mismatched points under RANSAC framework; 8 Compute its displacement using (3); 9 Estimate the optimal positions of all images by solving (4); 10 return the global panorama; For each pixel x i , it will be kept as one candidate crack pixel if the following criteria is satisfied, where g o (x i ) is the intensity of the original pixel, g t (x i ) is the intensity of corresponding pixel in the smoothed image, and T g is a specific threshold. Taking into account the different texture and illumination, the contrast between crack pixels and its surrounding background varies significantly, so the threshold T g is adaptively determined. To decease the effect of imhomogenous texture and illumination, we divide the whole image into grid cells, denoted as C i , i = 1, . . . , n c . In each grid cell C i , we select the threshold T g as where µ i and σ i are the mean and standard deviation of the difference between original intensities and their corresponding smoothed intensities of all pixels within C i , k g is an adjustment parameter which we generally set to be 1 in our application.
Since the image of airport runway may quite noisy, as illustrated in Fig. 3. The extracted candidate crack pixels based on intensity property contains many isolated or non-crack ones, which need to be found and removed according to their shape characteristic, as presented in the following step.

2) SHAPE BASED CRACK PIXEL FILTERING
The shape characteristic is used to distinguish the cracks from other dark patches in the image. First, we classify all candidate crack pixels into different groups, denoted as G i , i = 1, . . . , n g , according to their connectivity. To remove the noisy smaller blobs, we only keep the groups which contain more than T b pixels and whose contour/area ratio is large enough.
where notation | · | means the size of a collection, c(G i ) and a(G i ) denote the contour length and area of G i , respectively, T r is a specific threshold. While by applying the above process most crack pixels are detected, noises are also brought in, specifically linear pavement textures for concrete pavement, which possess the similar property, as shown in Fig. 3. However, pavement textures always appear as parallel straight lines, which makes them different from cracks. Thus, firstly, isolated pavement texture is detected as follows. At one hand, a minimum bounding box is detected for each candidate crack pixel group G i , that has a narrow width can be inferred as containing a straight line, which can be classified as pavement texture. At the other hand, pavement texture that is connected to cracks is detected as follows. Skeleton extraction based on K3M algorithm [27] is applied for each connected component G i to extract center axis of the patch, after which line detection using LSD [28] line segment detector is applied, and parallel straight lines can be classified as pavement textures. Morphological processing is applied afterwards to the detected pavement textures, taking into account that the pavement textures are always several pixels wide.
Leaving out the pavement texture pixels from the candidate crack pixels, we obtain the skeleton of cracks. Then we perform crack region growing by absorbing dark pixels neighboring the currently detected cracks, if the intensities of these pixels can satisfy a threshold test. The intensity threshold T n is determined as follows, where µ c and σ c are the mean and standard deviation of the grey levels of all currently detected crack pixels, and the parameter k c is set to be 1 in our experiments. This aggregation process is iteratively performed, so that all potential crack pixels surrounding the crack skeleton can be incorporated. Our proposed crack detection method is summarized in Algorithm 2. VOLUME 8, 2020

Algorithm 2 Crack Detection Algorithm
input : Image I output: all crack pixels in I 1 Divide I into grid cells C i , i = 1, . . . , n c ; 2 foreach C i do 3 Compute the intensity threshold T g using (6);

9
Remove G i if it belongs to a group of parallel line segments; 10 Extract the skeleton of G i using K3M algorithm; 11 Iteratively perform crack growing by absorbing dark pixels neighboring to G i if their intensities satisfy (8); 12 returnall crack pixels;

C. SUBSURFACE DEFECT DETECTION FROM GPR
As we discussed in section II, one of the main reasons, that existing CNN-based method cannot be directly deployed in applications, is that only 2D B-scan images are employed, which ignores the natural 3D property of subsurface defects. Thus, we propose an algorithm to handle this problem. Before introducing our algorithm, we first give some definitions and notations. While the GPR scans on the ground along a linear trajectory, it generates a B-scan image, that records the reflected pulses and shows the subsurface situation. When the GPR fixed on robot moves over a regular grid to collect multiple parallel B-scans, it produces a series of B-scans. This ensemble of B-scans forms a 3D data which is named as a C-scan. Denote B k as the k-th B-scan image, C = {B k |k = 1, . . . , n b } as the C-scan consisting of n b B-scans. We can evenly divide the 3D C-scan into slices of 2D images from main view, top view, and left view, respectively. The 2D images from main view are B-scan images. We denote T m , m = 1, . . . , n t and L n , n = 1, . . . , n l as the images sliced from top view and side view, respectively.
We propose a voting-based strategy to fuse the information of 2D images sliced from all views of C. We first adopt a CNN-based 2D object detection algorithm, such as Faster R-CNN, to generate the 2D bounding boxes labeled with categories and probabilities as the candidate defect regions. Denote r j as the j-th bounding box, with its probability belonging to category c k as p j,k . For each point X i in C, we define its score belonging to category c k to be the sum of probabilities from all 2D bounding boxes of B k , T m and L n . Thus, we compute the score as with where X i ∈ r j indicates point X i being inside of r j .
With the scores of all points in C computed, we only keep the points if the following criteria is satisfied.
where T s is a specific threshold for point scores.
Then, we cluster all remained points according to their positions and residing categories. DBSCAN algorithm [29] is adopted for this clustering step. DBSCAN starts with a remained point X i of category c k , and retrieves all points of the same category density-reachable from X i with regard to two parameters, ε and MinPts. The ε-neighborhood of X i is defined as, (12) where N ε (X i ) contains at least MinPts of remained points. In our experiments, we set ε = 16 and MinPts = 4.
Finally, for each point cluster, we generate a 3D bounding box as the detected defect region.
We summarize our subsurface defect detection method in Algorithm 3.

Algorithm 3 Subsurface Defect Detection Algorithm
input : C output: 3D bounding boxes of defects 1 Divide C into a set of B k , T m and L n ; 2 Generate candidate 2D bounding boxes from B k , T m and L n using Faster R-CNN; 3 foreach X i do 4 Compute its score s(X i , c k ) using (9); 5 Filter all points in C using (11); 6 Cluster the remain points according to their spacial intensities and categories using (12); 7 Generate 3D bounding box for each point cluster; 8 return all 3D bounding boxes;

V. EXPERIMENTS
We have employed our ARIR system in over 20 airports of China. To evaluate our proposed algorithms thoroughly, we select a representative airport of southeast of China for both quantitative and qualitative analysis. The surveyed region is over 2000 m 2 , and subsurface depth is 1.53 meters. Two human experts have labelled the surface cracks and subsurface defects individually. With the collected inspection data, we fulfill to perform both the surface and subsurface defect detection and visualization.

A. LARGE SCALE IMAGE STITCHING RESULTS
We succeed to stitch the large amount of images of the airport runways. An example of image stitching results of a region of 20m × 405m, with a total of 11642 original images is shown in Fig. 4, where we can see that the panorama of the surveyed airport runway is generated.

B. SURFACE CRACK DETECTION RESULTS
To quantitatively evaluate the performance of our crack detection method, we select the representative images to form a dataset, consisting of 112 images with the resolution of 1200 × 900 pixels. We employ three well-known metrics, including precision, recall and F1-measure, for the evaluation. The ground truth is manually labelled by human carefully. Considering it is inevitable to introduce errors when labelling, we assume that true positive pixels are included within a 5-pixel vicinity of the ground truth. These three metrics can be computed based on true positive (TP), false negative (FN), and false positive (FP), We present some example crack detection results in Fig. 5, where we can see that the images are quite challenging due to the special texture, the thin width of cracks and poor lighting conditions. The statistical quantitatively results are shown in Tab. 1. The results of our proposed algorithm can satisfy the requirement of field applications.

C. SUBSURFACE DEFECT DETECTION RESULTS
We compare our algorithm with Faster R-CNN which is only conducted on B-scan. To evaluate the performance of different methods quantitatively, three metrics, including Precision, Recall and F1-measure, are employed. To illustrate these metrics for object-level detection, we employ the indicator IoU which represents the overlapping rate of overlap between the resulting box and the labelled box. If the IoU value is greater than a specific threshold T IoU , we consider this box as a true positive. The quantitative results are shown in Tab. 2, where we can see that our proposed algorithm achieves a better performance compared with standard Faster R-CNN. Fig. 6 presents an example of defect qualitative results from two different views of a GPR C-scan segment. This global subsurface defect map can be used for quantitative analysis and also for visually assessing global defect patterns.

VI. CONCLUSION
We developed the automated airport runway inspection data collection and analysis. Specifically, we proposed a crack detection method to detect the surface cracks robustly from the collected images. A deep learning algorithm for GPR data analysis was presented to perform the subsurface defect detection. The image stitching algorithm allowed to generate an acceptable global image of the large scale airport runway. A data analysis software embedding these algorithms was developed. We tested our proposed system in applications over 20 airports and achieved satisfying performance.