Tunnel lining detection image void area identification method based on ground penetrating radar

Abstract

The invention discloses a tunnel lining detection image void area identification method based on a ground penetrating radar. The method comprises the following steps: collecting image data of a tunnel lining based on the ground penetrating radar; preprocessing the image data through a 1D-2D combined filtering processing method; carrying out image gray level enhancement processing on the preprocessed image based on gray level morphology and mathematical morphology operation; based on a multi-detection algorithm, edge detection is carried out on the image after gray level enhancement processing; carrying out adaptive threshold segmentation on the image after edge detection based on a maximum between-class variance method; and based on a connected region marking algorithm, removing small connected regions in the image after adaptive threshold segmentation. According to the recognition method, void area recognition can be automatically performed on a large number of ground penetrating radar images, the automation degree is high, and the influence of human factors on the recognition result of the void area in the tunnel lining can be effectively reduced. The method is simple in principle, the void area in the ground penetrating radar detection image can be identified more accurately, and the detection precision and efficiency are improved.

Description

Technical Field

[0001] The invention relates to the technical field of tunnel lining image recognition, and in particular to a method for identifying a void area in a tunnel lining detection image based on a ground penetrating radar. Background Art

[0002] Tunnel engineering is an important infrastructure construction project in the fields of transportation, water conservancy engineering, etc. During the construction and use of tunnels, the integrity and safety of the lining are crucial to the normal operation of the tunnel. However, due to factors such as construction quality, geological conditions, and service life, quality problems such as voids, cracks, and leakage may occur in the tunnel lining. These problems will not only affect the normal use of the tunnel, but may also pose a threat to personnel safety.

[0003] A major problem in the prior art is that the image features of the void detection behind the tunnel lining are not significant, which leads to low accuracy in void positioning and range identification. On the radar image, the void area may be difficult to distinguish from the surrounding rock mass, resulting in difficulty in identification. In addition, the void in the lining appears as a bright reflection different from the surrounding rock mass on the ground penetrating radar image, with an obvious arc-shaped reflection interface and messy internal reflections. Its image features are easily confused with pre-buried pipes. The current interpretation method of radar images is mostly manual experience interpretation of data, which has the disadvantages of low efficiency, strong subjectivity, and high requirements for personnel. This means that the detection results rely heavily on the experience and skills of the operator, lacking objectivity and consistency.

[0004] Ground penetrating radar is a non-destructive testing technology based on the principle of electromagnetic wave reflection. It is widely used in the detection of underground structures, rocks, soil and other geological bodies. It has also become a commonly used non-destructive testing method for tunnel lining detection. In order to solve the problem of hollow areas in tunnel lining, ground penetrating radar is used to conduct regular inspections of tunnel linings, timely identify hollow areas and repair problems, which is of great significance to ensure the safe operation of tunnels. Summary of the invention

[0005] The present invention provides a method for identifying hollow areas in tunnel lining inspection images based on ground penetrating radar. The method is based on the principle of image grayscale distribution law and performs hollow area identification through steps such as image acquisition, preprocessing, grayscale distribution analysis, and threshold setting. Hollow area identification can be automatically performed on a large number of ground penetrating radar images, with a high degree of automation and the ability to effectively reduce the influence of human factors on the identification results of hollow areas in tunnel linings. The method is simple in principle, can more accurately identify hollow areas in ground penetrating radar inspection images, and improve the accuracy and efficiency of detection.

[0006] The technical solution adopted by the present invention is:

[0007] The method for identifying the void area of ​​a tunnel lining inspection image based on ground penetrating radar includes the following steps:

[0008] Step 1: Collect image data of tunnel lining based on ground penetrating radar;

[0009] Step 2: Preprocess the image data by 1D-2D joint filtering method;

[0010] Step 3: Perform image grayscale enhancement processing on the image preprocessed in step 2 based on grayscale morphology and mathematical morphology operations;

[0011] Step 4: Perform edge detection on the image processed by grayscale enhancement in step 3 based on a multi-detection algorithm;

[0012] Step 5: Adaptive threshold segmentation of the image after edge detection in step 4 based on the maximum inter-class variance method;

[0013] Step 6: Based on the connected component labeling algorithm, remove the small connected domains in the image after the adaptive threshold segmentation in step 5.

[0014] In step 1, the ground penetrating radar (GPR) mainly uses the high-frequency pulse electromagnetic waves emitted by the antenna to detect the target. The CO1760 integrated ground penetrating radar is used as a detection device to continuously measure the tunnel lining structure through the test wheel and time trigger mode, and the point trigger mode is used for intermittent measurement. The test theory mode is used to detect the flat surface. When the test surface is relatively rough, the time trigger mode is used. A certain time interval is set according to the moving speed of the equipment to automatically collect data. Image data such as Figure 8 shown.

[0015] The step 2 comprises the following steps:

[0016] S2.1: First, use one-dimensional filtering technology to process the image, use static correction method to identify and mark the direct wave or initial reflected wave in the image data, determine the time delay of the direct wave by analyzing the initial part of the image data, remove the time delay to eliminate the influence of the direct wave on the image data, and obtain the zero offset section. The radar signal S'(t) after static correction is expressed as:

[0017] S′(t)=S(tt d ) (1);

[0018] In formula (1), t represents the test time, t d represents the time delay of the direct wave, S(tt d ) represents the radar signal after the time delay is eliminated.

[0019] S2.2: Then, the de-tributary drift method is used to eliminate the low-frequency DC component (also called DC drift or baseline offset) caused by factors such as instrument stability, environmental changes, or underground medium characteristics, including:

[0020] First calculate the average value of the radar signal average value It represents the DC component in the radar signal, that is, the low-frequency, stable or slowly changing part of the signal; then, the DC component is subtracted from each value of the original radar signal to return the baseline of the signal to zero or a predetermined reference level, which is mathematically expressed as:

[0021]

[0022] In the above formula, S(n) represents the discrete signal data; N represents the total number of sampling points in the discrete signal; n represents the number of sampling points, and S′(n) represents the discrete signal data after removing DC drift; Represents the average value of discrete radar signals.

[0023] S2.3: Next, energy attenuation processing is performed, by analyzing the energy attenuation characteristics of the radar signal over time, using exponential and power rate functions based on signal amplitude to compensate for energy attenuation. The processed signal has a more uniform energy distribution.

[0024]

[0025] In formula (4), α and β represent attenuation coefficients; G0 represents the signal amplitude at the initial moment (or reference frequency), G(f) represents the signal amplitude at frequency f, and f represents the frequency; f β Indicates the decay rate of signal amplitude as the frequency changes.

[0026] S2.4: Use two-dimensional filtering technology to process the image, identify the background noise in the radar image through background removal technology, and use filters to remove low-frequency noise and irrelevant background signals from the radar image. The details are as follows:

[0027] Use the moving average filter to remove background noise, improve the signal-to-noise ratio of the image, and highlight the structural features. The mathematical expression of the moving average filter is:

[0028]

[0029] In formula (5), f′(x, y) is the pixel value of the post-image at the coordinate (x, y); f(x+i, y+j) is the pixel value of the original image at the coordinate (x+i, y+j); M is the size of the filter window; i and j are the displacements relative to the center of the window (x, y), and their range is from arrive

[0030] S2.5: Use bandpass filtering to remove high-frequency noise and low-frequency drift from the radar signal while retaining the signal within a specific frequency range, as follows:

[0031] Butterworth bandpass filtering is used to distinguish the frequency range of noise signals from the frequency range of valid signals, and determine the signal frequency range that needs to be retained, that is, the bandpass frequency. The transfer function of the Butterworth filter allows signals within a specific frequency range to pass while suppressing signals of other frequencies. The transfer function is expressed as follows:

[0032]

[0033] In formula (6), H(f) represents the transfer function value of the filter at frequency f; f represents the frequency; f refers to the cutoff frequency; q is the filter order;

[0034] S2.6: Use the sliding average method to smooth the radar image and reduce the random noise in the image. By calculating the average value in the neighborhood of each point in the image, the random fluctuation of the image is reduced, making the continuous structural features in the image smoother and more continuous. The details are as follows:

[0035] Select a sliding window to cover a certain number of sampling points in the radar image. Average the signals in the window to calculate a signal. Slide the window along the signal and repeatedly calculate the average value to generate a smooth signal. Sliding average can reduce random noise and make continuous structural features smoother. The sliding average filter is specifically expressed as follows:

[0036]

[0037] In formula (7): g(i) represents the gray value of the filtered image at the i-th pixel; It is the pixel value of the original image at the corresponding position in the sliding window.

[0038] In step 3, the bright details and dark details in the image are enhanced by top-hat transformation in grayscale morphological operation. First, the original image is opened, and then the result of the opening operation is subtracted from the original image to highlight those areas that are brighter than the surrounding background. This area is called white top.

[0039] Then perform a closing operation on the image and subtract the original image from the result of the closing operation to highlight those small areas that are darker than the surrounding areas. This area is called the black top hat.

[0040] Exponential transformation and power law transformation based on mathematical morphology theory can also enhance the grayscale of images and highlight the detailed features of the empty area. Three transformations are used to enhance the image, and then the processing effects of the three image enhancements are comprehensively compared from the perspectives of background grayscale, empty features, and empty area grayscale information, and the optimal image processing result is selected for the next step of image edge detection.

[0041] The step 3 comprises the following steps:

[0042] S3.1: Perform top-hat transformation on the preprocessed image to enhance the local contrast of the image and highlight the details. First, use the white top-hat transformation to highlight the bright details, and perform an opening operation on the original image. The opening operation is a process of first erosion and then expansion. Use the structural element B to operate on the image f. The mathematical expression of the erosion operation is:

[0043]

[0044] in: represents the result of the erosion operation; f represents the image, B represents the structural element, and D is the domain of the image. For each pixel x in the image, the erosion operation finds all y, y represents the pixel mapping obtained by the structural element operation of x, so that yx falls within the structural element B. The maximum value of f(y); f(y) represents the gray value of the original image at y. The expansion operation is to expand the erosion result outward and restore some of the eroded boundaries. The mathematical expression of the expansion operation is as follows:

[0045]

[0046] in: represents the result of the dilation operation, which is to find the minimum value of g(y) for all pixels y in the image so that yx falls within -B (i.e., the reflection of the structural element B). g(y) represents the value of the structural element B at position y. Subtract the result of the opening operation from the original image to get the white top hat, highlighting the area that is brighter than the surrounding area:

[0047]

[0048] in: It means that the expansion operation is performed on the result of the corrosion operation, and White Top-hat represents the white top-hat area.

[0049] Then use the black top hat transform to highlight the dark details, and perform a closing operation on the original image. The closing operation is a process of dilation followed by erosion, and the structural element B is also used:

[0050]

[0051] in: Represents the result of the dilation operation, a represents the pixel point, and the dilation operation is to find the minimum value of g(b) for all pixels b in the image so that ba falls within -B (that is, the reflection of the structure element B) through the inverse structure element operation; g(b) represents the value of the structure element B at position b.

[0052] Then erode the dilation result:

[0053]

[0054] in: Represents the result of the erosion operation. For each pixel a in the image, find all b such that ba falls within the structure element B, and the maximum value of f(b); f(b) represents the grayscale value of the original image at b.

[0055] The result of the closing operation is then subtracted from the original image to create a black top hat, highlighting areas that are darker than the surrounding area:

[0056]

[0057] in: represents the result of the dilation operation; It means that the result of the dilation operation is further corroded, and White Top-hat indicates the black top-hat area.

[0058] S3.2: Use the exponential transformation and power rate transformation method to enhance the gray value of the hollow feature area and weaken the gray value of the noise gray area, thereby increasing the contrast between the hollow area and the noise area. The exponential transformation is used to enhance the high gray area in the image while retaining the low gray area. For each pixel (x, y) in the image, an exponential function is applied according to its gray value I(x, y). The form of the exponential transformation is:

[0059] I′(x,y)=u·e v·I(x,y) (14);

[0060] Among them: I′(x,y) represents the grayscale value after transformation; u and v are parameters that control the degree of transformation; u determines the maximum grayscale value of the image after transformation to prevent overflow; v determines the rate of grayscale value growth.

[0061] The power law transform enhances the de-focusing features of radar images by adjusting the image contrast while reducing the impact of noise. For each pixel in the image, a power law function is applied to adjust its grayscale value. The form of the power law transform is:

[0062] I′(x,y)=c·I(x,y) γ (15);

[0063] Among them, I′(x,y) is the gray value after transformation, c is the normalization coefficient, which ensures that the gray value of the transformed image is within a reasonable range; γ is the index for controlling contrast;

[0064] S3.3: Comprehensively compare the processing effects of the three image enhancements of S3.1 and S3.2 from the perspectives of background grayscale, void features, and void area grayscale information, and select the optimal image processing result for the next step of image edge detection.

[0065] In step 4, after the radar detection image is enhanced in step 3, the pixel brightness information of the image is greatly changed, and a clear boundary is formed between the image and the surrounding area. Extracting the effective boundary of the image can better reflect the curve characteristics of the empty area. Four different edge detection algorithms, namely, Roberts crossover operator, Prewitt edge detection operator, Sobel edge detection operator and Canny edge detection algorithm, are used to perform edge detection on the target image respectively, and then the results of the four detection algorithms are compared and analyzed, and the result with the best comprehensive effect is selected for the next step of adaptive threshold segmentation.

[0066] 4.1: First, the Roberts crossover operator is used to detect the edge of the target image. The edge is found by calculating the difference on the pixel diagonal. A 2x2 template is used and the operation is performed in two directions (diagonal direction). The algorithm definition and template are as follows:

[0067] For each pixel f(x,y) in the image, the difference in two directions is calculated using the following template:

[0068]

[0069] Roberts x Represents the edge detection template for the horizontal direction.

[0070]

[0071] Roberts y Represents the edge detection template for the vertical direction.

[0072] For the horizontal direction, the differential G x The calculation of (x,y) is as follows:

[0073] G x (x,y)=f(x+1,y)-f(x-1,y) (18);

[0074] Among them, f(x+1,y) represents the pixel value of the image at (x+1,y); f(x-1,y) represents the pixel value of the image at (x-1,y);

[0075] For the vertical direction, the differential G y The calculation of (x,y) is as follows:

[0076] G y (x,y)=f(x,y+1)-f(x,y-1) (19);

[0077] Among them, f(x,y+1) represents the pixel value of the image at (x,y+1); f(x,y-1) represents the pixel value of the image at (x,y-1).

[0078] value

[0079] Take the maximum value of the difference in two directions as the approximation of the edge strength G(x,y):

[0080] G(x,y)=max(|G x (x,y)|,|G y (x,y)|) (20);

[0081] Among them, G x (x, y) represents the gradient component of the Roberts crossover operator at the position (x, y) along the x-axis direction; G y (x,y) represents the gradient component along the y-axis at the position (x,y);

[0082] 4..2: Use the Prewitt edge detection operator to detect the edge of the target image. The Prewitt operator adds a mean filter step based on the Roberts operator. It uses the grayscale difference between adjacent pixels above and below and adjacent pixels below and below to extract the vertical and horizontal boundaries of the image, and expands the vertical template into three rows and the horizontal template into three columns. The following is the algorithm

[0083] Definitions and templates:

[0084] Horizontal edge detection template Prewitt x It is expressed as follows:

[0085]

[0086] Vertical edge detection template Prewitt y It is expressed as follows:

[0087]

[0088] Transverse edge strength G x The calculation is as follows:

[0089] G x=[f(x-1,y-1)+f(x-1,y)+f(x-1,y+1)]-[f(x+1,y-1)+f(x+1,y)+f(x+1,y+1)](23);

[0090] Among them, f(x-1,y-1) represents the pixel value of the image at (x-1,y-1), f(x-1,y) represents the pixel value of the image at (x-1,y), f(x-1,y+1) represents the pixel value of the image at (x-1,y+1), f(x+1,y-1) represents the pixel value of the image at (x+1,y-1), f(x+1,y) represents the pixel value of the image at (x+1,y), and f(x+1,y+1) represents the pixel value of the image at (x+1,y+1);

[0091] Longitudinal edge strength G y The calculation is as follows:

[0092] G y =[f(x-1,y+1)+f(x,y+1)+f(x+1,y+1)]-[f(x-1,y-1)+f(x,y-1)+f(x+1,y-1)](24);

[0093] The total edge strength G(x,y) is calculated as follows

[0094]

[0095] Among them, G x (x,y) represents Prewitt x The gradient component of the crossover operator along the x-axis at the position (x, y); G y (x,y) represents the gradient component along the y-axis at the position (x,y);

[0096] 4.3: Use the Sobel edge detection operator to detect the edge of the target image, and calculate the edge strength by weighted summing the gray values ​​of the pixels. The algorithm definition and template are as follows:

[0097] Horizontal edge detection template Sobel x It is expressed as follows:

[0098]

[0099] Vertical edge detection template Sobel y It is expressed as follows:

[0100]

[0101] Transverse edge strength G x The calculation is as follows:

[0102] G x =[f(x+1,y-1)-f(x-1,y-1)]+2[f(x+1,y)-f(x-1,y)]+[f(x+1,y+1)-f(x-1,y-1)](28);

[0103] Longitudinal edge strength G y The calculation is as follows:

[0104] G y =[f(x-1,y+1)-f(x-1,y-1)]+2[f(x,y+1)-f(x,,y-1)]+[f(x+1,y+1)-f(x-1,y-1)] (29);

[0105] The total edge strength G(x,y) is calculated as follows:

[0106]

[0107] 4.4: Use the Canny edge detection algorithm to detect the edge of the target image, and use a variant of the Sobel operator to calculate the gradient magnitude G and direction of the image. The gradient magnitude and gradient direction are calculated using the following formula:

[0108]

[0109] Where: G x and G y are the gradient components in the horizontal and vertical directions respectively; G x and G y Calculate using the Sobel operator method in 4..3.

[0110] Then at each pixel position, check the gradient direction and take the maximum gradient value only in that direction. Set the high threshold T h and low threshold T l , for each pixel, according to its gradient magnitude G i For classification, if G i >T h , then the pixel is considered a strong edge pixel. i ≤T l , then the pixel is considered a non-edge pixel. l <G i ≤T h , then mark these pixels as weak edge pixels. Then initialize all pixels of the image to 0, set all strong edge pixels to 1 and mark them as edge pixels. For each weak edge pixel, if it is adjacent to a marked strong edge pixel, set its pixel to 1 and mark it as an edge pixel to form a continuous edge.

[0111] The image edge detection results obtained by the four detection methods are comprehensively compared and analyzed from the perspective of edge feature information and edge continuity, and the optimal image processing result is selected for the next step of image threshold segmentation.

[0112] In step 5, image segmentation is used to quickly screen the feature areas and invalid areas in the image and then extract the required target image. However, the data contained in the target area of ​​the ground penetrating radar image is too large, and a threshold needs to be set for image segmentation. In the actual detection process, the background noise of the ground penetrating radar image is large, and setting a single threshold cannot meet the requirements of displaying the image missing features. Therefore, adaptive threshold segmentation is used to compare the gray value difference of each pixel in the image.

[0113] The threshold is automatically obtained based on the maximum inter-class variance method. This method divides the image into two parts, foreground and background, according to the grayscale characteristics of the image. The foreground is the target feature. When the inter-class variance between the two parts is the largest, it means that the difference between the foreground and the background is the largest. The grayscale value at this time is the optimal threshold. After obtaining the optimal threshold, the image is segmented based on it. The grayscale value of each pixel in the image is compared with the threshold, and the pixel is set to black (0) or white (255) according to the comparison result; for each pixel in the grayscale image, if its grayscale value is less than or equal to the threshold, the pixel is set to 0 (black); indicating the background, if the grayscale value is greater than the threshold, the pixel is set to 255 (white), indicating the foreground. The details are as follows:

[0114] The threshold is automatically obtained based on the maximum inter-class variance method. First, the grayscale histogram H(i) of the image is calculated, where i is the grayscale level and H(i) is the number of times grayscale level i appears in the image.

[0115] Then the histogram is converted into a normalized density function P(i), which represents the probability of gray level i appearing in the entire image. P(i) is calculated as follows:

[0116]

[0117] Then, for each possible threshold t, the image is divided into foreground (grey value greater than or equal to t) and background (grey value less than t). The foreground is the target feature. The mean μ0 and μ1 of these two categories and their probabilities P0 and P1 are calculated as follows:

[0118]

[0119]

[0120] Among them, P0(t) and P1(t) are the probabilities of background and foreground under threshold t, respectively, and μ0(t) and μ1(t) are the corresponding means.

[0121] Then calculate the inter-class variance σ corresponding to the threshold t 2 (t):

[0122] σ 2 (t)=P0(t)·P1(t)·(μ0(t)-μ1(t)) 2 (38);

[0123] Then traverse all possible thresholds t and find the inter-class variance σ 2 (t) The largest t, denoted as t max :

[0124] t max = argmax t σ 2 (t) (39);

[0125] After obtaining the optimal threshold, image segmentation is performed based on it. The grayscale value of each pixel in the image is compared with the threshold, and the pixel is set to black (0) or white (255) according to the comparison result. For each pixel in the grayscale image, if its grayscale value is less than or equal to the threshold, the pixel is set to 0 (black); indicating the background, if the grayscale value is greater than the threshold, the pixel is set to 255 (white), indicating the foreground, as shown below:

[0126]

[0127] Among them, f(x,y) is the grayscale value of the original image at position (x,y), and E(x,y) is the pixel value of the edge image at position (x,y).

[0128] In step 6, the small connected domain in the image is a region composed of pixels with the same pixel value and adjacent positions in the image, that is, noise, details and unnecessary small regions in the image.

[0129] The connected component labeling algorithm is used to identify independent sets of pixels in the image that have the same pixel value g(x, y) and are adjacent to each other (i.e., connected domains).

[0130] After the connected regions are marked, the size of each connected domain (i.e., the number of pixels in it) is obtained. Then, by setting a threshold T, connected domains smaller than the threshold are regarded as "small connected domains" and deleted from the image (set to background pixel values).

[0131] After deleting the small connected domain, the noise and small, unimportant details in the image will be eliminated, while the larger, important objects or areas will be retained, and the final empty area image will be obtained, which is specifically expressed as follows:

[0132]

[0133] Among them, g′(x,y) represents the image pixel value after removing the small connected domain.

[0134] The present invention provides a method for identifying a void area in a tunnel lining inspection image based on a ground penetrating radar, and the technical effects are as follows:

[0135] 1) Advantages of step 1 of the present invention: Radar data of tunnel linings are collected based on ground penetrating radar non-destructive testing technology, which can detect tunnel linings without destroying the structure, detect deep structural problems and provide high-resolution images.

[0136] 2) Advantages of step 2 of the present invention: Image preprocessing of radar data based on the 1D-2D joint filtering processing method can effectively remove the noise introduced during the data acquisition process, improve the image quality, and provide clear data for subsequent image enhancement and edge detection.

[0137] 3) Advantages of step 3 of the present invention: grayscale enhancement processing of the image based on grayscale morphology and mathematical morphology operations can effectively highlight the detailed features of the cracks and void areas in the image.

[0138] 4) Advantages of step 4 of the present invention: Image edge detection based on multiple detection algorithms can effectively provide edge information at different angles, improve the robustness of detection and reduce possible misjudgments of a single algorithm.

[0139] 5) Advantages of step 5 of the present invention: Adaptive threshold segmentation of images based on the maximum inter-class variance method can automatically determine the threshold to reduce the intervention of human factors, improve consistency and repeatability, and thus obtain the best segmentation effect.

[0140] 6) Advantages of step 6 of the present invention: Removing small connected domains in the image based on the connected region labeling algorithm can effectively eliminate false features generated by noise or image processing, making the main features of the missing area more prominent. BRIEF DESCRIPTION OF THE DRAWINGS

[0141] The present invention will be further described below in conjunction with the accompanying drawings and examples;

[0142] Figure 1 This is the plan view of the physical model of the tunnel lining structure.

[0143] Figure 2 Comparison chart before and after image preprocessing.

[0144] Figure 3 This is the image enhancement effect diagram.

[0145] Figure 4 This is the image edge detection result diagram.

[0146] Figure 5 This is the threshold segmentation effect diagram.

[0147] Figure 6 Comparison diagram before and after removing small connected domains.

[0148] Figure 7 Extract images for the voided areas.

[0149] Figure 8 Image of the tunnel lining as detected by ground penetrating radar. DETAILED DESCRIPTION

[0150] The method for identifying the void area of ​​a ground penetrating radar tunnel lining inspection image based on grayscale distribution includes the following steps:

[0151] Step 1: Collect radar data of tunnel lining based on ground penetrating radar non-destructive testing technology;

[0152] Step 2: Image preprocessing of radar data based on 1D-2D joint filtering method;

[0153] Step 3: Perform image grayscale enhancement processing based on grayscale morphology and mathematical morphology operations;

[0154] Step 4: Perform image edge detection based on multi-detection algorithm;

[0155] Step 5: Perform image adaptive threshold segmentation based on the maximum inter-class variance method;

[0156] Step 6: Remove small connected domains in the image based on the connected component labeling algorithm.

[0157] In the step 1, radar data of the tunnel lining is collected based on the ground penetrating radar nondestructive testing technology, including the following steps:

[0158] The CO1760 integrated ground penetrating radar is used as a detection device to continuously measure the physical model of the constructed tunnel lining structure. The detection target depth is 0.5m-1m, and the detection medium is engineering concrete. The ground penetrating radar is parameterized and the distance measuring wheel mode is used to make the radar close to the model surface at a certain speed to detect at a constant speed, so as to obtain a stable distance and time depth profile.

[0159] In the step 2, image preprocessing is performed on the radar data based on the 1D-2D joint filtering method, including the following steps:

[0160] Step S2.1: Use static correction method to identify and mark the direct wave or initial reflected wave in the radar signal. Determine the time delay t of the direct wave by analyzing the initial part of the original radar signal S(t). d, subtract this time delay from the radar signal to eliminate the influence of the direct wave, thereby obtaining a zero offset profile. The statically corrected S'(t) is expressed as:

[0161] S′(t)=S(tt d )

[0162] Step S2.2: Use the DC drift removal method to eliminate the low-frequency DC component (also called DC drift or baseline offset) caused by factors such as instrument stability, environmental changes, or underground medium characteristics. First, calculate the average value of the radar signal The resulting mean represents the DC component of the signal, that is, the low-frequency, stable or slowly changing part of the signal. The DC component is then subtracted from each value of the original signal to return the baseline of the signal to zero or a predetermined reference level. The mathematical expression of this step is:

[0163]

[0164] Where N is the total number of samples in the signal.

[0165] Step S2.3: Analyze the energy attenuation characteristics of the radar signal over time. Apply a gain function G(f), usually an exponential or power-rate function based on the signal amplitude, to compensate for the energy attenuation. The processed signal has a more uniform energy distribution.

[0166] G(f)=e -αf (Exponential decay)

[0167]

[0168] Where α and β represent the attenuation coefficients, and f represents the frequency.

[0169] Step S2.4: Use a moving average filter to remove background noise, improve the signal-to-noise ratio of the image, and highlight the structural features. The mathematical expression of the moving average filter is:

[0170]

[0171] Where M is the window size of the filter.

[0172] Step S2.5: Use Butterworth bandpass filtering to distinguish the frequency range of the noise signal and the frequency range of the valid signal, and determine the signal frequency range that needs to be retained, that is, the bandpass frequency. The transfer function of the Butterworth filter allows signals within a specific frequency range to pass, while suppressing signals of other frequencies. Its transfer function is expressed as follows:

[0173]

[0174] Where f is the frequency, f is the cutoff frequency, and N is the filter order.

[0175] Step S2.6: Select a sliding window to cover a certain number of sampling points in the radar image. Average the signals within the window to calculate a signal. Slide the window along the signal and repeatedly calculate the average to generate a smooth signal. Sliding average can reduce random noise and make continuous structural features smoother.

[0176]

[0177] In the step 3, the image grayscale enhancement processing is performed based on grayscale morphology and mathematical morphology operations, including the following steps:

[0178] Step S3.1: Perform top-hat transformation on the preprocessed image to enhance the local contrast of the image and highlight the details. First, use the white top-hat transformation to highlight the bright details, and perform an opening operation on the original image. The opening operation is a process of first erosion and then expansion. Use the structural element B to operate on the image f. The mathematical expression of the erosion operation is:

[0179]

[0180] Where D is the domain of the image, and B is the translation of the structural element B centered on x. The dilation operation is to expand the eroded result outward and restore some of the eroded boundaries:

[0181]

[0182] Subtracting the result of the opening operation from the original image yields a white top hat that highlights areas that are brighter than their surroundings:

[0183]

[0184] Then use the black top hat transform to highlight the dark details, and perform a closing operation on the original image. The closing operation is a process of dilation followed by erosion, and the structural element B is also used:

[0185]

[0186] Then erode the dilation result:

[0187]

[0188] Subtracting the result of the closing operation from the original image produces a black top hat that highlights areas that are darker than their surroundings:

[0189]

[0190] Step S3.2: Use the exponential transformation and power rate transformation method to enhance the gray value of the empty feature area and weaken the gray value of the noise gray area, thereby increasing the contrast between the empty area and the noise area. The exponential transformation is used to enhance the high gray area in the image while retaining the low gray area. For each pixel in the image, an exponential function is applied according to its gray value I(x,y). The form of the exponential transformation is:

[0191] I′(x,y)=a·e b·I(x,y)

[0192] I′(x,y) represents the grayscale value after transformation, and a and b are parameters that control the degree of transformation. a determines the maximum grayscale value of the image after transformation to prevent overflow; b determines the rate at which the grayscale value increases.

[0193] The power law transform enhances the de-focusing features of radar images by adjusting the image contrast while reducing the impact of noise. For each pixel in the image, a power law function is applied to adjust its grayscale value. The form of the power law transform is:

[0194] I′(x,y)=c·I(x,y) γ

[0195] Among them, I′(x, y) is the gray value after transformation, c is the normalization coefficient, which ensures that the gray value of the transformed image is in a reasonable range, and γ is the index for controlling the contrast.

[0196] Step S3.3: From the processing effects of the three image enhancements in steps 3.1 and 3.2, the top hat transformation operator evenly processes the background grayscale, and the contrast of the missing features is not obvious. The exponential transformation and power rate transformation have better enhancement effects, and the grayscale information brightness of the missing area of ​​the exponential transformation is high, which can meet the requirements of further processing. The image processing result of the exponential transformation is selected for the next step of image edge detection.

[0197] In the step 4, image edge detection is performed based on a multi-detection algorithm, including the following steps:

[0198] Step S4.1: After the radar detection image is enhanced in step 3, the pixel brightness information of the image has changed greatly, and a clear boundary is formed between the image and the surrounding area. Extracting the effective boundary of the image can better reflect the curve characteristics of the empty area. First, the Roberts crossover operator is used to detect the edge of the target image. The edge is found by calculating the difference on the pixel diagonal. A 2x2 template is used and the operation is performed in two directions (diagonal direction). The algorithm definition and template are as follows:

[0199] G(x,y)=|f(x,y)-f(x+1,y+1)|+|f(x,y+1)-f(x+1,y)|

[0200]

[0201] Step S4.2: Use the Prewitt edge detection operator to perform edge detection on the target image. The Prewitt operator adds a mean filter step based on the Roberts operator, and uses the grayscale difference between adjacent pixels above and below and adjacent pixels below and below to extract the vertical and horizontal boundaries of the image, expanding the vertical template into three rows and the horizontal template into three columns. The following is the algorithm definition and template:

[0202] G x =[f(x-1,y-1)+f(x-1,y)+f(x-1,y+1)]-[f(x+1,y-1)+f(x+1,y)+f(x+1,y+1)]

[0203] G y =[f(x-1,y+1)+f(x,y+1)+f(x+1,y+1)]-[f(x-1,y-1)+f(x,y-1)+f(x+1,y-1)]

[0204]

[0205] Step S4.3: Use the Sobel edge detection operator to perform edge detection on the target image, and calculate the edge strength by weighted summing the gray values ​​of the pixels. The algorithm definition and template are as follows:

[0206] d x =[f(x+1,y-1)-f(x-1,y-1)]+2[f(x+1,y)-f(x-1,y)]+[f(x+1,y+1)-f(x-1,y-1)]

[0207] d y =[f(x-1,y+1)-f(x-1,y-1)]+2[f(x,y+1)-f(x,,y-1)]+[f(x+1,y+1)-f(x-1,y-1)]

[0208]

[0209] Step S4.4: Use the Canny edge detection algorithm to perform edge detection on the target image, and use a variant of the Sobel operator to calculate the gradient magnitude G and direction of the image. The gradient magnitude is calculated using the following formula:

[0210]

[0211] Among them G x and G yare the gradient components in the horizontal and vertical directions respectively.

[0212] Then at each pixel position, check the gradient direction and take the maximum gradient value only in that direction. Set a high threshold and a low threshold respectively. The high threshold is used to determine strong edges, and pixels below the low threshold are marked as weak edges. Finally, these edge pixels are tracked to form a continuous edge.

[0213] Step S4.5: Comprehensively compare and analyze the image edge detection results obtained by the four detection methods from the perspective of edge feature information and edge continuity. All four detection methods can effectively remove background noise, but compared with the three processing results of Roberts, Prewitt and Canny operators, Sobel's processing result has better boundary characteristics, so the Sobel operator detection result is selected for further threshold segmentation.

[0214] In the step 5, image adaptive threshold segmentation is performed based on the maximum inter-class variance method, including the following steps:

[0215] Based on the maximum inter-class variance method, the threshold is automatically obtained. According to the grayscale characteristics of the image, the image is divided into two parts: foreground and background. The foreground is the target feature. When the inter-class variance between the two parts is the largest, it means that the difference between the foreground and the background is the largest. At this time, the grayscale value is the optimal threshold. After obtaining the optimal threshold, the image is segmented according to it. The grayscale value of each pixel in the image is compared with the threshold, and the pixel is set to black (0) or white (255) according to the comparison result. For each pixel in the grayscale image, if its grayscale value is less than or equal to the threshold, the pixel is set to 0 (black); indicating the background, if the grayscale value is greater than the threshold, the pixel is set to 255 (white), indicating the foreground.

[0216] In step 6, removing small connected domains in the image based on a connected domain labeling algorithm includes the following steps:

[0217] The connected domain labeling algorithm is used to identify independent sets of pixels in the image that have the same pixel value and are adjacent to each other (i.e., connected domains). After the connected domains are labeled, the size of each connected domain (i.e., the number of pixels in it) is obtained. Then, by setting a threshold, connected domains smaller than the threshold are regarded as "small connected domains" and deleted from the image (usually set to the background pixel value). After deleting small connected domains, noise and small, unimportant details in the image will be eliminated, while larger, important objects or areas will be retained to obtain the final empty area image.

[0218] Figure 2The figure is a comparison of the image before and after preprocessing. Since the DC drift removal, background removal, bandpass filtering and sliding average processing in the preprocessing can reduce the random noise in the image and improve the signal-to-noise ratio, the preprocessed image shows less noise than before preprocessing, and the preprocessed image is clearer than the original image because the filtering process helps to highlight the structural features and make the lining structure in the image more obvious. The energy gain step compensates for the energy attenuation of the radar signal, making the preprocessed image more uniform in energy distribution and improving the consistency of the signal.

[0219] Figure 3 The image enhancement effect is shown in Figure 2. It can be seen that the exponential transformation enhances the high grayscale areas in the image while retaining the low grayscale areas. The power rate transformation enhances the hollow features by adjusting the image contrast while reducing the noise impact. These two transformations help improve the visibility of the hollow areas in the image, making the subsequent edge detection and feature recognition more accurate.

[0220] Figure 4 is the image edge detection result. Figure 4 It can be seen that the Roberts operator uses the difference in the diagonal direction to find the edge, showing a more obvious diagonal edge feature in the image. The Prewitt operator extracts the edge by calculating the grayscale difference between the upper and lower adjacent pixels and the left and right adjacent pixels, showing a smoother edge contour in the image. The Sobel operator calculates the edge strength by weighted summation of the grayscale values ​​of the pixels, providing clearer edge information in the image, especially in the strength and direction of the edge. The Canny operator is a multi-stage algorithm that uses high and low thresholds to determine strong and weak edges, showing continuous and less noisy edges in the image.

[0221] Figure 5 This is the effect diagram of threshold segmentation. Figure 5 It can be seen that the empty area (foreground) in the image is separated from the background through threshold segmentation, and the contrast between the features of the empty area and the background is stronger, making it easier to be identified.

[0222] Figure 6 This is a comparison chart before and after removing small connected domains. Figure 6 It can be seen that after removing the small connected domain, the main features (empty areas) in the image are clearer because irrelevant noise and small details have been eliminated.

[0223] Figure 7 Extract images for the empty areas. Figure 7 The voids in the tunnel lining are clearly visible and highlighted in the image for easy identification.

Claims

1. A method for identifying void areas in tunnel lining inspection images based on ground penetrating radar, characterized in that The following steps are involved: Step 1: Collect image data of tunnel lining based on ground penetrating radar; Step 2: Preprocess the image data by 1D-2D joint filtering method; Step 3: Perform image grayscale enhancement processing on the image preprocessed in step 2 based on grayscale morphology and mathematical morphology operations; Step 4: Perform edge detection on the image processed by grayscale enhancement in step 3 based on the detection algorithm; Step 5: Adaptive threshold segmentation of the image after edge detection in step 4 based on the maximum inter-class variance method; Step 6: Based on the connected component labeling algorithm, remove the small connected domains in the image after the adaptive threshold segmentation in step 5.

2. The method for identifying void areas in tunnel lining inspection images based on ground penetrating radar according to claim 1 is characterized by: The step 2 comprises the following steps: S2.1: First, the static correction method is used to identify and mark the direct wave or initial reflected wave in the image data. The time delay of the direct wave is determined by analyzing the initial part of the image data. The time delay is removed to eliminate the influence of the direct wave on the image data, thereby obtaining the zero offset section. The radar signal S'(t) after static correction is expressed as: S′(t)=S(t-t d )(1); In formula (1), t represents the test time, t d represents the time delay of the direct wave, S(tt d ) represents the radar signal after the time delay has been eliminated; S2.2: Then, the low-frequency DC component is eliminated by using the branch drift removal method, which includes: First calculate the average value of the radar signal average value It represents the DC component in the radar signal, that is, the low-frequency, stable or slowly changing part of the signal; then, the DC component is subtracted from each value of the original radar signal to return the baseline of the signal to zero or a predetermined reference level, which is mathematically expressed as: In the above formula, S(n) represents the discrete signal data; N represents the total number of sampling points in the discrete signal; n represents the number of sampling points, and S′(n) represents the discrete signal data after removing DC drift; represents the average value of discrete radar signals; S2.3: Next, energy attenuation processing is performed, by analyzing the energy attenuation characteristics of the radar signal over time, and using exponential and power rate functions based on the signal amplitude to compensate for the energy attenuation; In formula (4), α and β represent attenuation coefficients; G0 represents the signal amplitude at the initial moment (or reference frequency), G(f) represents the signal amplitude at frequency f, and f represents the frequency; f β Indicates the attenuation rate of signal amplitude as the frequency changes; S2.4: Use a moving average filter to remove background noise, improve the signal-to-noise ratio of the image, and highlight the structural features. The mathematical expression of the moving average filter is: In formula (5), f′(x, y) is the pixel value of the post-image at the coordinate (x, y); f(x+i, y+j) is the pixel value of the original image at the coordinate (x+i, y+j); M is the size of the filter window; i and j are the displacements relative to the center of the window (x, y), and their range is from arrive S2.5: Use Butterworth bandpass filtering to distinguish the frequency range of noise signals and the frequency range of valid signals, and determine the signal frequency range that needs to be retained, that is, the bandpass frequency; the transfer function of the Butterworth filter allows signals within a specific frequency range to pass, while suppressing signals of other frequencies; the transfer function is expressed as follows: In formula (6), H(f) represents the transfer function value of the filter at frequency f; f represents the frequency; f refers to the cutoff frequency; q is the filter order; S2.6: Select a sliding window to cover a certain number of sampling points in the radar image; average the signals in the window to calculate a signal; slide the window along the signal and repeatedly calculate the average value to generate a smooth signal; the sliding average filter is specifically expressed as follows: In formula (7): g(i) represents the gray value of the filtered image at the i-th pixel; It is the pixel value of the original image at the corresponding position in the sliding window.

3. The method for identifying void areas in tunnel lining inspection images based on ground penetrating radar according to claim 1 is characterized by: The step 3 comprises the following steps: S3.1: Perform top-hat transformation on the preprocessed image to enhance the local contrast of the image and highlight the details. First, use white top-hat transformation to highlight the bright details, and perform an opening operation on the original image. The opening operation is a process of first erosion and then expansion. Use the structure element B to operate on the image f. The mathematical expression of the erosion operation is: in: represents the result of the erosion operation; f represents the image, B represents the structural element, and D is the domain of the image. For each pixel x in the image, the erosion operation finds all y, y represents the pixel mapping obtained by the structural element operation of x, so that yx falls within the structural element B. The maximum value of f(y); f(y) represents the gray value of the original image at y; the expansion operation is to expand the erosion result outward and restore some of the eroded boundaries. The mathematical expression of the expansion operation is as follows: in: represents the result of the dilation operation. The dilation operation is to find the minimum value of g(y) for all pixel points y in the image so that yx falls within -B. g(y) represents the value of the structural element B at position y. Subtracting the result of the opening operation from the original image yields a white top hat that highlights areas that are brighter than their surroundings: in: Indicates that the result of the corrosion operation is then expanded. White Top-hat indicates the white top-hat area. Then use the black top hat transform to highlight the dark details, and perform a closing operation on the original image. The closing operation is a process of dilation followed by erosion, and the structural element B is also used: in: represents the result of the dilation operation, a represents the pixel point, and the dilation operation is to find the minimum value of g(b) for all pixels b in the image so that ba falls within -B through the inverse structure element operation; g(b) represents the value of the structure element B at position b; Then erode the dilation result: in: represents the result of the erosion operation. For each pixel a in the image, find all b such that ba falls within the structure element B, and the maximum value of f(b); f(b) represents the gray value of the original image at b. The result of the closing operation is then subtracted from the original image to create a black top hat, highlighting areas that are darker than the surrounding area: in: represents the result of the dilation operation; Indicates that the result of the dilation operation is subjected to the erosion operation, and White Top-hat indicates the black top-hat area; S3.2: The exponential transformation and power rate transformation methods are used to enhance the grayscale value of the hollow feature area and weaken the grayscale value of the noise grayscale area, thereby increasing the contrast between the hollow area and the noise area; the exponential transformation is used to enhance the high grayscale area in the image while retaining the low grayscale area; for each pixel point (x, y) in the image, the exponential function is applied according to its grayscale value I(x, y), and the form of the exponential transformation is: I′(x,y)=u·e v·I(x,y) (14); Where: I′(x,y) represents the gray value after transformation; u and v are parameters that control the degree of transformation; u determines the maximum gray value of the image after transformation to prevent overflow; v determines the rate of gray value growth; The power law transform enhances the clearance features of radar images by adjusting the image contrast while reducing the influence of noise. For each pixel in the image, a power law function is applied to adjust its gray value. The form of the power law transform is: I′(x,y)=c·I(x,y) γ (15); Among them, I′(x,y) is the gray value after transformation, c is the normalization coefficient, which ensures that the gray value of the transformed image is within a reasonable range; γ is the index for controlling contrast; S3.3: Comprehensively compare the processing effects of the three image enhancements of S3.1 and S3.2 from the perspectives of background grayscale, void features, and void area grayscale information, and select the optimal image processing result for the next step of image edge detection.

4. The method for identifying void areas in tunnel lining inspection images based on ground penetrating radar according to claim 1 is characterized by: In step 4, the Sobel edge detection operator is used to perform edge detection on the target image, and the edge strength is calculated by weighted summing of the gray values ​​of the pixels. The algorithm definition and template are as follows: Horizontal edge detection template Sobel x It is expressed as follows: Vertical edge detection template Sobel y It is expressed as follows: Transverse edge strength G x The calculation is as follows: G x =[f(x+1,y-1)-f(x-1,y-1)]+2[f(x+1,y)-f(x-1,y)]+[f(x+1,y+1)-f(x-1,y-1)](28); Longitudinal edge strength G y The calculation is as follows: G y =[f(x-1,y+1)-f(x-1,y-1)]+2[f(x,y+1)-f(x,,y-1)]+[f(x+1,y+1)-f(x-1,y-1)](29); The total edge strength G(x,y) is calculated as follows:

5. The method for identifying void areas in tunnel lining inspection images based on ground penetrating radar according to claim 1 is characterized by: In step 5, the threshold is automatically obtained based on the maximum inter-class variance method, and the grayscale histogram H(i) of the image is first calculated, where i is the grayscale level and H(i) is the number of times the grayscale level i appears in the image; Then the histogram is converted into a normalized density function P(i), which represents the probability of gray level i appearing in the entire image. P(i) is calculated as follows: Then, for each possible threshold t, the image is divided into foreground and background. The foreground is the target feature. Calculate the means μ0 and μ1 of the two categories and their probabilities P0 and P1 as follows: Among them, P0(t) and P1(t) are the probabilities of background and foreground under threshold t, respectively, and μ0(t) and μ1(t) are the corresponding means; Then calculate the inter-class variance σ corresponding to the threshold t 2 (t): s 2 (t)=P0(t)·P1(t)·(μ0(t)-μ1(t)) 2 (38); Then traverse all possible thresholds t and find the inter-class variance σ 2 (t) The largest t, denoted as t max : t max =argmax t σ 2 (t)(39); After obtaining the optimal threshold, image segmentation is performed based on it. The grayscale value of each pixel in the image is compared with the threshold, and the pixel is set to black (0) or white (255) according to the comparison result. For each pixel in the grayscale image, if its grayscale value is less than or equal to the threshold, the pixel is set to 0 (black); indicating the background, if the grayscale value is greater than the threshold, the pixel is set to 255 (white), indicating the foreground, as shown below: Among them, f(x,y) is the grayscale value of the original image at position (x,y), and E(x,y) is the pixel value of the edge image at position (x,y).

6. The method for identifying void areas in tunnel lining inspection images based on ground penetrating radar according to claim 1, characterized in that: In step 6, a connected region labeling algorithm is used to identify a set of pixels in the image that are independent of each other, have the same pixel value g(x, y) and are adjacent in position, i.e., a connected region; After the connected regions are marked, the size of each connected domain, i.e., the number of pixels in it, is obtained; then, by setting a threshold T, connected domains smaller than the threshold are regarded as "small connected domains" and deleted from the image; After deleting the small connected domain, the noise and small, unimportant details in the image will be eliminated, while the larger, important objects or areas will be retained, and the final empty area image will be obtained, which is specifically expressed as follows: Among them, g′(x,y) represents the image pixel value after removing the small connected domain.

Images