A method and system for the joint detection of internal and surface defects in a road using ground penetrating radar

Abstract

The application discloses a kind of inside and surface defect joint method and system of fusing ground penetrating radar for road, belong to road detection field, by collecting surface defect and internal radar data and unified mileage coordinate;Extract surface crack center, width, strike and radar anomaly body center, depth, reflection intensity;By electromagnetic wave propagation and solid mechanics coupling forward to generate simulation data set, establish to the nonlinear mapping of defect dip angle and horizontal offset;The measured characteristics are input into the mapping, and the predicted dip angle and offset are output, if the offset is greater than the set proportion of depth, it is judged as non-vertical correlation, otherwise it is judged as vertical correlation;Finally, mark the correlation type in three-dimensional visualization;The application can accurately identify inclined defects, eliminate homing error, and avoid false grouting and collapse accidents.

Description

Technical Field

[0001] This invention relates to the field of road inspection, and more specifically to a method and system for combining ground-penetrating radar to detect internal and surface defects in roads. Background Technology

[0002] In existing technologies, the joint detection of internal road defects and surface cracks usually involves using ground-penetrating radar to obtain information on underground anomalies, and then combining this with surface images to identify cracks. The judgment is made using the "vertical nearest correlation" criterion (that is, the radar anomaly directly below the surface crack is considered to be the cause of the crack).

[0003] However, this method does not take into account complex geological scenarios such as the development of tilted defects or horizontal offsets. Furthermore, conventional offset imaging tends to misplace the reflected energy of tilted interfaces directly below cracks, resulting in misplacement errors. This can lead to misclassifying non-vertical correlations as vertical correlations, which in turn can cause erroneous grouting that avoids the real defects or even road collapse accidents, seriously affecting the accuracy of maintenance decisions and the safety of road operation. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for integrating ground-penetrating radar to detect road interior and surface defects, thereby addressing the shortcomings of the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for combining ground-penetrating radar with analysis of road interior and surface defects, comprising: Simultaneously collect road surface defect data and internal radar data along the road detection direction; Each frame of surface defect data and each internal radar data is assigned a unified mileage coordinate, so that the same mileage position corresponds to both a surface defect data unit and an internal radar data unit. Image segmentation and edge detection are performed on the surface defect data unit to extract the center coordinates, width, and direction of the surface cracks; offset imaging and anomaly identification are performed on the internal radar data unit to extract the center coordinates, burial depth, and reflection intensity of the radar anomaly. A simulation dataset with multiple combinations of defect parameters is generated in advance through the coupled forward modeling of electromagnetic wave propagation and solid mechanical response. Based on the simulation dataset, a nonlinear mapping function is established from the center coordinates, burial depth, reflection intensity of radar anomaly, and the center coordinates, width, and orientation of surface cracks to the defect inclination angle and horizontal offset. The center coordinates, width, and orientation of the extracted surface cracks, as well as the center coordinates, burial depth, and reflection intensity of the radar anomaly, are input into the nonlinear mapping function, which outputs the predicted defect tilt angle and the predicted horizontal offset. If the predicted horizontal offset is greater than a set ratio of the burial depth, it is determined to be a non-vertical correlation caused by radar offset repositioning error; otherwise, it is determined to be a vertical correlation. Based on the judgment results, the association type between internal road defects and surface cracks is generated, and three types of states—vertical association, non-vertical association, or uncertain—are marked in the road 3D visualization data.

[0006] Preferably, the surface defect data includes image information and elevation information of the road surface, and the internal radar data is electromagnetic wave reflection signal collected by a multi-channel or array antenna.

[0007] Preferably, each set of defect parameters includes the defect tilt angle and the horizontal offset between the defect and the surface crack, and each set of simulation datasets generates a radar anomaly coordinate and burial depth as well as a surface crack coordinate.

[0008] Preferably, the specific steps for performing image segmentation and edge detection on the surface defect data unit to extract the center coordinates, width, and orientation of the surface crack are as follows: Anisotropic diffusion filtering is performed on the grayscale image in the surface defect data unit to preserve crack edges and suppress road texture noise; The second-order partial derivative matrix of each pixel is calculated for the filtered image. Linear structure pixels are selected based on the ratio of the two eigenvalues ​​of the second-order partial derivative matrix to form crack candidate regions. Morphological skeleton extraction is performed on the candidate region to obtain the center line of the crack with a width of one pixel. The gray-level gradient is calculated along the normal direction of the center line, and the gradient peak spacing is used as the crack width. The centerline is segmented and fitted with straight lines, and the direction angle of the fitted straight lines is used as the crack direction.

[0009] Preferably, the specific steps for performing offset imaging and anomaly identification on the internal radar data unit to extract the center coordinates, burial depth, and reflection intensity of the radar anomaly are as follows: Time zero-point correction and exponential gain compensation are performed on the eight waveform sequences in each internal radar data unit. The energy of each reflected wave is reversed along the propagation path and superimposed to generate a two-dimensional depth profile with depth as the vertical axis and lateral distance as the horizontal axis. The ratio of local energy to background noise energy is calculated point by point on the depth profile, and continuous pixel regions with a ratio greater than 4 are marked as candidate anomalous bodies. Calculate the lateral coordinates and depth values ​​of the geometric center of the candidate anomaly on the depth profile, and use them as the center coordinates and burial depth. Take the maximum amplitude of all pixels in the region as the reflection intensity.

[0010] Preferably, the forward modeling of the coupling between electromagnetic wave propagation and solid mechanical response includes: The first step is to construct a three-dimensional mesh array that includes dielectric constant and elastic modulus assignments, and set the defect tilt angle and horizontal offset in the mesh; The second step involves using the finite-difference time-domain method to solve the electromagnetic wave propagation equation on the same grid array to obtain the radar waveform sequence; then, the finite element method is used to solve the stress balance equation to obtain the surface displacement field. The third step is to store the defect tilt angle, horizontal offset, anomaly coordinates and burial depth in the radar waveform sequence, and crack coordinates and width in the surface displacement field as a simulated data record according to their corresponding relationships. The fourth step involves iterating through the defect inclination angle from 0 degrees to 60 degrees and the horizontal offset from 0 meters to 1.5 meters, repeating steps two and three to generate a simulation dataset covering all parameter combinations.

[0011] Preferably, establishing the nonlinear mapping function and outputting the predicted defect tilt angle and predicted horizontal offset includes: The input parameters for each set of data are extracted from the simulated dataset, including: radar anomaly center coordinates, burial depth, reflection intensity, surface crack center coordinates, width, and orientation; and the output parameters include: defect tilt angle and horizontal offset. Using Gaussian radial basis functions, with the input parameter vector of each sample point as the center, the Euclidean distance between any two sample point input vectors is calculated, and a system of linear equations of order with the number of sample points is constructed. Solve the system of linear equations to obtain the radial basis function weight coefficients for each sample point; For the measured input parameters, calculate the Euclidean distance between them and the input vector of each sample point, substitute them into the radial basis function, and then sum them with weighted coefficients to obtain the predicted defect tilt angle and the predicted horizontal offset.

[0012] Preferably, for each detection section, the association type at the section is generated based on the judgment result: the judgment result is obtained by comparing the predicted horizontal offset with the burial depth: if the predicted horizontal offset is greater than 0.3 times the burial depth, the association type is marked as non-vertical association; if the predicted horizontal offset is less than or equal to 0.3 times the burial depth, the association type is marked as vertical association; if the confidence level of the radial basis function interpolation output is less than 0.7 and a judgment cannot be made, the association type is marked as uncertain.

[0013] Preferably, three-dimensional visualization data of the road is constructed, with the starting and ending mileage of the detected road segment as the longitudinal range, the road cross-section width as the lateral range, and the detection depth as the vertical range, to establish a three-dimensional spatial grid; each grid point is assigned an initial color value, where grid points above the road surface are set to gray, and grid points below the road surface without anomalies are set to semi-transparent blue.

[0014] This invention also provides a combined system for integrating ground-penetrating radar with road interior and surface defects, comprising: Data acquisition module: Simultaneously acquires road surface defect data and internal radar data along the road detection direction; Spatial registration module: Assigns unified mileage coordinates to each frame of surface defect data and each internal radar data, so that the same mileage position corresponds to one surface defect data unit and one internal radar data unit at the same time. Feature extraction module: performs image segmentation and edge detection on the surface defect data unit to extract the center coordinates, width and direction of surface cracks; performs offset imaging and anomaly recognition on the internal radar data unit to extract the center coordinates, burial depth and reflection intensity of radar anomalies; Forward modeling and mapping module: In advance, through the coupled forward modeling of electromagnetic wave propagation and solid mechanical response, a simulation dataset under multiple combinations of defect parameters is generated, and based on the simulation dataset, a nonlinear mapping function is established from the center coordinates, burial depth, reflection intensity of the radar anomaly, and the center coordinates, width, and orientation of the surface crack to the defect inclination angle and horizontal offset. The correlation judgment module inputs the center coordinates, width, and direction of the extracted surface cracks, as well as the center coordinates, burial depth, and reflection intensity of the radar anomaly into the nonlinear mapping function, and outputs the predicted defect tilt angle and the predicted horizontal offset. If the predicted horizontal offset is greater than the set ratio of the burial depth, it is determined to be a non-vertical correlation caused by the radar offset positioning error; otherwise, it is determined to be a vertical correlation. Visualization output module: Based on the judgment results, it generates the association type between internal road defects and surface cracks, and marks three types of states in the road 3D visualization data: vertical association, non-vertical association, or uncertain.

[0015] The technical effects and advantages provided by the present invention in the above technical solution are as follows: 1. This invention overcomes the limitations of existing "vertical nearest correlation" criteria in terms of physical mechanism. Existing technologies only assume that defects develop vertically upwards, ignoring the asymmetry of radar wave propagation paths in scenarios with tilted defects or horizontal offsets. This leads to conventional offset imaging incorrectly repositioning the reflected energy of the tilted interface directly below the crack, resulting in repositioning errors. This invention, through coupled forward modeling of electromagnetic wave propagation and solid mechanical response, simultaneously simulates the electromagnetic wave reflection response to the tilted defect and the mechanical deformation response of the road surface above the defect under load in the same three-dimensional mesh array. This reveals how the defect tilt angle and horizontal offset respectively affect the lateral repositioning position of the radar anomaly and the spatial distribution of surface cracks. Based on this coupled physical relationship, a nonlinear mapping function can directly inversely deduce the true tilt angle and horizontal offset of the defect from the measured radar anomaly coordinates, burial depth, reflection intensity, and surface crack coordinates, width, and orientation, eliminating the interference of offset repositioning errors on correlation judgment from a physical perspective.

[0016] 2. This invention directly solves the technical problem of incorrect grouting and road collapse caused by misjudgment. Traditional methods misjudge inclined defects with large horizontal offsets as vertically related, leading maintenance personnel to grout directly below the crack (where there is actually no defect). This not only fails to fill the true defect but may also accelerate collapse due to the disturbance caused by the grouting pressure. This invention, by predicting the ratio of horizontal offset to burial depth, clearly determines whether it is a non-vertical correlation, thus indicating the actual spatial location of the true defect. This determination result is directly translated into a red sphere (the true location of the anomaly) and a red dashed line (offset relationship) in the visual annotation, enabling maintenance personnel to accurately perform grouting at the true center of the defect. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0018] Figure 1 This is a flowchart of a method for combining ground-penetrating radar with methods for detecting internal and surface defects in roads, according to the present invention.

[0019] Figure 2 This is a flowchart of a combined system module for integrating ground-penetrating radar and road interior and surface defects according to the present invention.

[0020] Figure 3 This is a flowchart of the forward modeling method for the coupling of electromagnetic wave propagation and solid mechanical response of the present invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] Example 1, please refer to Figure 1 As shown in this embodiment, a method for combining ground-penetrating radar to analyze road interior and surface defects includes: Road surface defect data and internal radar data are collected simultaneously along the road detection direction. The surface defect data includes image information and elevation information of the road surface, and the internal radar data consists of electromagnetic wave reflection signals collected by multi-channel or array antennas.

[0023] First, the surface data acquisition component and the internal radar data acquisition component need to be installed on the same inspection vehicle.

[0024] The surface data acquisition component employs a linear array image capture device and a laser elevation measurement device. The linear array image capture device has at least 2048 pixel units, and its capture frequency is set to trigger once every 0.01 meters of travel to acquire a grayscale image of the road surface. The laser elevation measurement device uses the principle of multi-point laser ranging, with 16 laser ranging probes evenly arranged laterally. Each probe also records the vertical distance every 0.01 meters, thereby obtaining the elevation values ​​of the road surface at 0.05-meter intervals laterally.

[0025] The internal radar data acquisition component uses a linear array consisting of eight transceiver radar antenna elements. The antenna elements are arranged at equal intervals in the horizontal direction, with a spacing of 0.1 meters between adjacent antenna elements. The operating center frequency is 900 MHz. Each antenna element can simultaneously transmit and receive electromagnetic wave reflection signals.

[0026] The detection vehicle travels at a constant speed of 30 kilometers per hour. To ensure strict synchronization between surface data and radar data in terms of mileage, a rotating coded ranging wheel is installed at the rear wheel of the vehicle. Each rotation of this ranging wheel outputs 2000 pulse signals, corresponding to 0.1 meters of vehicle movement. Each pulse signal simultaneously triggers the linear array image acquisition device, the laser elevation measurement device, and the radar antenna array to acquire data. Thus, for every 0.1 meters of movement, a set of surface grayscale image data, a set of elevation data at 16 points laterally, and a set of electromagnetic wave reflection waveform data recorded by each of the eight radar antenna elements are obtained.

[0027] For surface defect data, grayscale images acquired by the linear array imaging device are transmitted to the storage unit via data cable and continuously saved in the format of "mileage value - image data block". The 16 elevation values ​​acquired by the laser elevation measurement device are also stored in mileage order and linked with grayscale images of the same mileage. For internal radar data, the reflected waveforms acquired by each radar antenna element are recorded with 512 sampling points per channel and a time window length of 40 nanoseconds. Each sampling point value represents the instantaneous amplitude of the echo signal. Data from all eight antenna elements are stored in the storage unit in the format of "mileage value - channel number - waveform sampling sequence".

[0028] Each frame of surface defect data and each internal radar data is assigned a unified mileage coordinate, so that the same mileage location corresponds to both a surface defect data unit and an internal radar data unit.

[0029] Before the test begins, the pulse counter of the rotary encoder distance measuring wheel installed at the rear wheel of the test vehicle is reset to zero, and the actual road mileage marker value of the starting point is manually recorded. For example, if the starting point marker is "12 km + 345 m", the value in meters is 12345 meters. This starting point mileage value serves as the benchmark for all subsequent mileage coordinate calculations.

[0030] After the vehicle starts and travels at a constant speed, the ranging wheel outputs 2000 pulse signals per revolution, corresponding to the vehicle moving 0.1 meters forward. Each pulse signal simultaneously triggers the linear array image capture device, the laser elevation measurement device, and the radar antenna array to collect data, and records the current cumulative pulse count as the nth pulse, where n increments from 1. At this point, the mileage coordinates of the current data acquisition section are calculated based on the pulse number n: the mileage coordinates equal to the starting mileage value plus n multiplied by 0.1 meters.

[0031] For example, if the starting mileage is 12345 meters, when the 100th pulse arrives, the mileage coordinate is equal to 12345 plus 100 multiplied by 0.1, which is 12355 meters. This mileage coordinate is the unified spatial identifier for all surface defect data and internal radar data collected under this pulse trigger.

[0032] For surface defect data, when each pulse is triggered, the linear array image capturing device acquires a grayscale image with a width of 2048 pixels, and the laser elevation measuring device acquires a set of elevation values ​​at 16 lateral positions.

[0033] The grayscale image frame, the set of 16 elevation values, and the currently calculated mileage coordinates are stored together as a surface defect data unit. The specific storage format is: one record contains a mileage coordinate field, a grayscale image data block field, and 16 floating-point elevation value fields. For internal radar data, each pulse triggers an independent electromagnetic wave reflection waveform with 512 sampling points per pulse, and each waveform contains 512 instantaneous amplitude values. All eight waveform data collected under the same pulse are stored together with the same mileage coordinate as an internal radar data unit. The storage format is: one record contains a mileage coordinate field, eight waveform data block fields, and each waveform data block is a sequence of 512 amplitude values.

[0034] During the detection process, a surface defect data unit and an internal radar data unit are generated for every 0.1 meters of advance, both with identical mileage coordinates. For example, at the 100th pulse, the mileage coordinates of both the surface defect data unit and the internal radar data unit are 12355 meters.

[0035] For any given mileage coordinate, the corresponding surface defect data unit and internal radar data unit can be uniquely indexed through that coordinate value, achieving precise spatial alignment between the two.

[0036] If abnormalities such as slippage of the measuring wheel or loss of pulse signals occur during the detection process, linear interpolation will be used to complete the mileage coordinates.

[0037] For example: Record the mileage coordinates and pulse number of the last valid pulse before the anomaly occurred, and the mileage coordinates and pulse number of the first valid pulse after the anomaly recovery. Assume the valid pulse number before the anomaly is 'a', and its mileage coordinate is 'La'; the valid pulse number after the anomaly is 'b', and its mileage coordinate is 'Lb'. For lost pulses with serial numbers between a+1 and b-1, their mileage coordinates are distributed at equal intervals, i.e., the mileage coordinate of the kth lost pulse is equal to La plus (k minus a) multiplied by (Lb minus La) divided by (b minus a). This calculation process ensures that the mileage coordinates at the lost section remain uniformly distributed, thus maintaining consistency with the coordinates of adjacent normal sections.

[0038] Using the above method, after the detection is completed, a set of continuous surface defect data unit sequences and a set of continuous internal radar data unit sequences are obtained. The two are in one-to-one correspondence according to the mileage coordinates, providing an accurate spatial reference for subsequent preprocessing and feature extraction.

[0039] Taking a real road section as an example: for a 500-meter-long section, a total of 5,000 pulse cross-sections were collected. The mileage coordinates of each cross-section range from the starting point of 12,345 meters to the ending point of 12,845 meters, with a step size of 0.1 meters. Under each mileage coordinate, there is simultaneously a surface defect data unit containing grayscale image and elevation data, as well as an internal radar data unit containing 8 radar waveforms.

[0040] Image segmentation and edge detection are performed on the surface defect data unit to extract the center coordinates, width, and direction of surface cracks; offset imaging and anomaly identification are performed on the internal radar data unit to extract the center coordinates, burial depth, and reflection intensity of radar anomalies.

[0041] In this embodiment of the invention, the specific implementation method for performing image segmentation and edge detection on the grayscale image in the surface defect data unit to extract the center coordinates, width and direction of the surface crack is as follows.

[0042] Anisotropic diffusion filtering is performed on the grayscale image. The core of this filtering method is that the grayscale value of each pixel is iteratively updated according to the diffusion coefficient of its grayscale difference with neighboring pixels. When the grayscale difference is large, the diffusion coefficient approaches 0, thus preserving crack edges; when the grayscale difference is small, the diffusion coefficient approaches 1, thus smoothing road surface texture noise. The number of iterations is set to 10. In each iteration, the updated value of each pixel is equal to its original value plus the sum of the weighted diffusion amounts in four directions (up, down, left, and right) of its neighboring pixels. The weights are determined by the absolute value of the grayscale difference in that direction. In the filtered image, crack edges are clear, while fine granular noise is suppressed.

[0043] Calculate the second-order partial derivative matrix for each pixel in the filtered image. The second-order partial derivative matrix is ​​also called the Hessian matrix. For each pixel in the image, calculate its second-order partial derivative along the horizontal direction, its second-order partial derivative along the vertical direction, and its mixed second-order partial derivative along both directions.

[0044] The specific calculation method is as follows: For the second-order partial derivative in the horizontal direction, the sum of the gray values ​​of the pixel's right-hand neighbor and left-hand neighbor is subtracted from twice the pixel's gray value, and then divided by the square of the distance between adjacent pixels. Similarly, the second-order partial derivative in the vertical direction is calculated using the pixels above and below each other. The mixed second-order partial derivative is calculated using the combination of pixels at the four corners: top left, bottom right, top right, and bottom left. After obtaining the three second-order partial derivative values ​​for each pixel, a 2x2 matrix is ​​formed. The two eigenvalues ​​of this matrix are calculated, and the ratio of the larger eigenvalue to the smaller eigenvalue is called the eigenvalue ratio. For linear structures (such as cracks), this ratio is much greater than 1; for speckled noise or uniform regions, this ratio is close to 1. A threshold of 5 is set for the ratio, and pixels with a ratio greater than 5 are marked as linear structure pixels; these pixels constitute the crack candidate region.

[0045] Morphological skeleton extraction is performed on the candidate region. Specifically, pixels at the edges of the candidate region are successively stripped, but those pixels whose removal would alter the region's connectivity are retained, until the remaining pixel width is a single pixel. The resulting single-pixel width curve is the crack centerline. Along the normal direction of each pixel on this centerline (perpendicular to the tangent direction of the centerline), a grayscale profile spanning the crack is extracted from the original grayscale image. On this profile, the grayscale values ​​at the edges of the crack show a sharp drop and rise. The grayscale difference between adjacent pixels on the profile is calculated, and the points of fastest grayscale decrease and rise are found. The pixel distance between these two points is the crack width at that location.

[0046] The average width of the crack is obtained by taking the arithmetic mean of the width values ​​of all pixels along the center line.

[0047] The center line of the crack is divided into several segments of equal length, each segment being approximately 20 pixels long.

[0048] For the centerline pixel coordinates within each segment, a straight line is fitted using the least squares method. The direction angle of this line (the angle rotated counterclockwise with the horizontal direction as 0 degrees) is the direction of the crack in that segment. If the direction of the entire crack does not change much, the average of the direction angles of all segments is taken as the overall direction.

[0049] In this embodiment of the invention, offset imaging and anomaly identification are performed on the 8-channel waveform sequence in the internal radar data unit to extract the center coordinates, burial depth, and reflection intensity of the radar anomaly. Specifically, this includes: Time zero-point correction and exponential gain compensation are performed on the eight waveform sequences in each internal radar data unit. The time zero-point correction method is as follows: the average amplitude of the first 10 sampling points of each waveform sequence is taken as the background noise level. The first sampling point in the waveform sequence that exceeds the background noise level by more than 5 times is determined as the arrival time of the ground reflection. The time coordinate of this sampling point is forcibly set to 0 nanoseconds, and the time coordinates of all subsequent sampling points are shifted forward in sequence.

[0050] The exponential gain compensation method works as follows: For the corrected waveform sequence, the amplitude value of each sampling point is multiplied by a coefficient that increases exponentially with the sampling point number. The base of this coefficient is a natural constant, and the exponent is the sampling point number multiplied by 0.01. After compensation, the amplitude of the weak signal in the deeper layers (within a large time window) is amplified to a level comparable to that of the shallower signals.

[0051] The energy of each reflected wave is reversed along the propagation path and superimposed to generate a two-dimensional depth profile with depth as the vertical axis and lateral distance as the horizontal axis.

[0052] Starting from the waveform sequence recorded at each receiving antenna location, based on the propagation speed of electromagnetic waves in the road material (set to 0.1 meters per nanosecond), the amplitude value of each sampling point is backpropagated to all possible reflection points underground. For a certain depth point underground, the algebraic sum of the backpropagation amplitudes of all antennas at that point is calculated as the offset imaging value for that point. The imaging values ​​are arranged according to depth (from 0 meters to 3 meters, with intervals of 0.01 meters) and lateral distance (from 0 meters to 0.7 meters, corresponding to the lateral range covered by 8 antennas, with intervals of 0.01 meters) to form a two-dimensional depth profile.

[0053] Calculate the ratio of local energy to background noise energy point-by-point on the depth profile. Local energy is defined as the sum of the squares of the amplitudes of all pixels within a circular region with a radius of 3 pixels centered at that point. Background noise energy is defined as the sum of the squares of the amplitudes of all pixels across the entire depth profile, divided by the total number of pixels, and then multiplied by the area of ​​the local region. Calculate the ratio of local energy to background noise energy for each pixel, and mark consecutive pixel regions with a ratio greater than 4 as candidate anomalous bodies.

[0054] For example, on a certain depth profile, if the ratio of eight consecutive pixels located at a horizontal distance of 0.35 meters and a depth of 1.2 meters is greater than 4, then these eight pixels constitute a candidate anomalous body.

[0055] Calculate the lateral coordinates and depth values ​​of the geometric center of the candidate anomaly on the depth profile, and use them as the center coordinates and burial depth.

[0056] The geometric center coordinates are calculated as follows: the arithmetic mean of the horizontal coordinates of all pixels within the candidate anomaly body is taken to obtain the center horizontal coordinates; the arithmetic mean of the depth values ​​of all pixels is taken to obtain the center depth value, which is the burial depth. The maximum amplitude of all pixels within the candidate anomaly body is taken as the reflection intensity of the anomaly.

[0057] Taking the above 8 candidate anomalies as an example, assuming their horizontal coordinates are 0.34 meters, 0.35 meters, 0.36 meters, etc., with an average of 0.35 meters; and their depth values ​​are 1.19 meters, 1.20 meters, 1.21 meters, etc., with an average of 1.20 meters, then the center coordinates of the anomaly are 0.35 meters horizontally and 1.20 meters deep, with a burial depth of 1.20 meters. The reflection intensity is taken as the maximum amplitude among these 8 pixels.

[0058] It should be noted that the grayscale image in the surface defect data unit effectively preserves the crack edges and suppresses road surface texture noise through anisotropic diffusion filtering. Then, linear structure pixels are selected based on the ratio of eigenvalues ​​of the second-order partial derivative matrix. The center line of a single pixel is obtained by combining morphological skeleton extraction. The crack width is obtained by calculating the gradient peak spacing along the normal direction. The direction is obtained by fitting the center line segmentally with straight lines, thus achieving accurate quantification of crack location, width and direction. At the same time, the 8 waveform sequences in the internal radar data unit are enhanced by time zero-point correction and exponential gain compensation. The reflected wave energy is reversed to generate a two-dimensional depth profile. Candidate anomalies are marked by the ratio of local energy to background noise energy (threshold 4). The geometric center coordinates and maximum amplitude are used as the center coordinates, burial depth and reflection intensity, thus effectively identifying the spatial location and signal strength of underground anomalies.

[0059] Please see Figure 3 As shown, a simulation dataset with multiple sets of defect parameters is generated in advance through the coupling forward modeling of electromagnetic wave propagation and solid mechanical response. Each set of defect parameters includes the defect tilt angle and the horizontal offset between the defect and the surface crack. Each set of simulation datasets corresponds to the coordinates and burial depth of a radar anomaly and the coordinates of a surface crack.

[0060] In this embodiment of the invention, the specific implementation method for generating a simulation dataset under multiple combinations of defect parameters in advance through the coupling forward modeling of electromagnetic wave propagation and solid mechanical response is as follows.

[0061] A three-dimensional grid array was constructed. The grid array has a horizontal length of 0.7 meters, consistent with the horizontal coverage of the eight radar antennas during detection; a horizontal length of 0.1 meters in the vertical direction (i.e., the direction of road movement), the same as the spacing between two adjacent data acquisitions; and a vertical length of 3.0 meters, covering the typical detection depth from the road surface to the roadbed. The grid spacing in all three directions was set to 0.01 meters, resulting in 71 horizontal grid points, 11 vertical grid points, and 301 depth grid points, for a total of 71 x 11 x 301 = 235,081 grid points.

[0062] Two physical properties are assigned to each grid point: dielectric constant (dimensionless, used for electromagnetic wave propagation calculations) and elastic modulus (in megapascals, used for solid mechanics calculations). The dielectric constant is assigned layer by layer according to the road structure: 6 for the surface layer (depth 0 to 0.2 meters), 8 for the base layer (depth 0.2 to 0.6 meters), and 12 for the subgrade (depth 0.6 to 3.0 meters). The elastic modulus is also assigned layer by layer in the same way: 3000 MPa for the surface layer, 800 MPa for the base layer, and 50 MPa for the subgrade. Defect regions are set in the grid array, with ellipsoidal shapes. Their inclination angle (the angle between the major axis of the defect and the horizontal plane) and horizontal offset (the horizontal distance between the defect center and the crack) are used as variable parameters.

[0063] The dielectric constant within the defect region is changed to 4 (representing the dielectric constant of air in the cavity), and the elastic modulus is changed to 0.1 MPa (representing that the cavity is nearly stiff). The defect tilt angle ranges from 0 degrees to 60 degrees, and the horizontal offset ranges from 0 meters to 1.5 meters.

[0064] On the same three-dimensional mesh array, the electromagnetic wave propagation equation is first solved using the finite-difference time-domain method. The core of the finite-difference time-domain method is to approximate the time and space partial derivatives in Maxwell's equations using finite difference forms.

[0065] The time axis is discretized into time steps, with a time step size of 0.01 nanoseconds. Within each time step, all grid points are traversed, and the electric and magnetic field values ​​for the next moment are calculated using the difference formula based on the electric and magnetic field values ​​of the current grid point and its adjacent grid points. A radar transmitting antenna is simulated on the upper surface of the grid array (at a depth of 0 meters), and a Ricker wavelet pulse with a center frequency of 900 MHz is applied as the excitation source. A receiving point is placed at a horizontal position of 0.35 meters on the upper surface of the grid array (i.e., the 4th of the 8 antennas), and the electric field intensity value at this point is recorded at each time step. Recording is performed continuously for 4000 time steps, corresponding to a total time window of 40 nanoseconds, resulting in an electromagnetic wave reflection waveform with a length of 4000 sampling points. The above process is repeated for the receiving point position, traversing 8 horizontal positions (0.1-meter intervals, from 0 meters to 0.7 meters), to obtain an 8-channel waveform sequence.

[0066] Identifying reflected waves from anomalies in waveform sequences: Find the peak point with the largest absolute amplitude in the waveform. Multiply the two-way travel time corresponding to this peak point by the propagation speed of electromagnetic waves in road materials (0.1 m / s) and divide by 2 to obtain the burial depth of the anomaly. The lateral coordinate of the receiving point corresponding to this peak point is the lateral coordinate of the center of the anomaly.

[0067] The stress balance equations are solved using the finite element method on the same three-dimensional mesh array.

[0068] Using a mesh array as the computational domain, a uniformly distributed load of 0.7 MPa (simulating a standard axle load) is applied to the upper surface. The bottom boundary of the mesh array is fixed (zero displacement), and the normal displacement of the four side boundaries is zero. Tetrahedral elements are used to mesh the surface, with the displacement field assumed to be linearly distributed within each element. The global stiffness matrix and load vector are assembled, and the linear equations are solved to obtain the displacement value of each mesh node. The vertical displacement components of each node on the upper surface are extracted to obtain the surface displacement field. The location with the largest change in vertical displacement gradient within the surface displacement field is identified; this location indicates the location where cracks will appear.

[0069] Calculate the vertical displacement difference between adjacent nodes on the upper surface, and take the line connecting the adjacent nodes with the largest difference as the crack location. The coordinates of the midpoint of this line are taken as the crack center coordinates. The crack width is obtained by converting the displacement difference with the material fracture toughness: it is set that when the vertical displacement difference exceeds 0.5 mm, the crack is considered to have opened, and the opening amount is the crack width.

[0070] The calculated defect dip angle, horizontal offset, transverse coordinates of the anomaly center, anomaly burial depth, crack center coordinates, and crack width are stored as a single simulation data record according to their corresponding relationships. The defect dip angle is iterated from 0 degrees to 60 degrees, with a value every 5 degrees; the horizontal offset is iterated from 0 meters to 1.5 meters, with a value every 0.1 meters. A total of 13 x 16 = 208 combinations are performed, with each combination repeating the electromagnetic wave and mechanical calculations, generating 208 simulation data records, forming a complete simulation dataset.

[0071] Each record contains the following fields: defect inclination angle (degrees), horizontal offset (meters), transverse coordinates of the anomaly center (meters), anomaly burial depth (meters), crack center coordinates (meters), and crack width (millimeters).

[0072] Based on the simulated dataset, a nonlinear mapping function is established from the center coordinates, burial depth, and reflection intensity of the radar anomaly, as well as the center coordinates, width, and orientation of the surface crack, to the defect tilt angle and horizontal offset. The extracted center coordinates, width, and orientation of the surface crack, along with the center coordinates, burial depth, and reflection intensity of the radar anomaly, are input into the nonlinear mapping function, which outputs the predicted defect tilt angle and the predicted horizontal offset. If the predicted horizontal offset is greater than a predetermined proportion of the burial depth, it is determined to be a non-vertical correlation caused by radar offset repositioning error; otherwise, it is determined to be a vertical correlation.

[0073] First, the input and output parameters for each set of data are extracted from the simulated dataset.

[0074] The simulation dataset contains 208 records, each corresponding to a set of defect parameter combinations. For each record, seven input parameters are extracted: the lateral coordinates of the radar anomaly center (in meters), the radar anomaly burial depth (in meters), the radar anomaly reflection intensity (dimensionless, ranging from 0 to 1), the lateral coordinates of the surface crack center (in meters), the surface crack width (in millimeters), and the surface crack orientation (in degrees, measured counterclockwise with the horizontal direction as 0 degrees). Two output parameters are also extracted: the defect tilt angle (in degrees, ranging from 0 to 60 degrees) and the defect horizontal offset (in meters, ranging from 0 to 1.5 meters). The input parameters from these 208 records are combined into a 208-row, 7-column input matrix, and the output parameters are combined into a 208-row, 2-column output matrix.

[0075] Gaussian radial basis functions are used to establish the mapping relationship between input and output parameters. The specific form of the Gaussian radial basis function is: ;in, Indicates sample points and σ represents the radial basis function value of the independent variable; exp represents the exponential function with the natural constant e≈2.71828 as the base; σ represents the width parameter of the radial basis function (also known as the smoothing factor); the width parameter is 0.5 times the maximum value of the Euclidean distance between all sample points.

[0076] Constructing a system of linear equations: Let the unknown radial basis function weight coefficient matrix be W, with a size of 208 rows and 2 columns. Then the system of linear equations can be expressed as: the radial basis function matrix R multiplied by the weight coefficient matrix W equals the output parameter matrix Y (208 rows and 2 columns). Here, the first column of the i-th row of the output parameter matrix Y represents the defect tilt angle of the i-th sample point, and the second column represents the horizontal offset of the i-th sample point.

[0077] Solving this system of linear equations yields the weight coefficient matrix W. Since the radial basis function matrix R is a symmetric positive definite matrix, Gaussian elimination is used to solve it.

[0078] The augmented matrix (R and Y combined) is transformed into an upper triangular matrix through row operations, and then each unknown is solved by back substitution. The solution process is performed column by column: first, the first column (the weight coefficients corresponding to the defect tilt angle) is solved, and then the second column (the weight coefficients corresponding to the horizontal offset) is solved. Two sets of weight coefficients are obtained corresponding to 208 sample points: the first set of 208 weight coefficients is used to predict the defect tilt angle, and the second set of 208 weight coefficients is used to predict the horizontal offset.

[0079] For the seven input parameters obtained from actual measurements, the same radial basis function interpolation method is used to predict the defect tilt angle and horizontal offset.

[0080] Let the measured input vector be... It includes the measured transverse coordinates of the radar anomaly center, burial depth, reflection intensity, transverse coordinates of the surface crack center, width, and orientation.

[0081] For the measured input vector Input vector of the i-th sample point Euclidean distance Defined as: ; Calculate the radial basis function values, radial basis function values Using Gaussian functions: ; Value range (0, 1).

[0082] Predicting defect tilt angle and horizontal offset, including: Predicted value of defect tilt angle for: ; Predicted value of horizontal offset for: ; in and These are the defect tilt angle weighting coefficient and the horizontal offset weighting coefficient, respectively, obtained by solving a system of linear equations, corresponding to the i-th sample point.

[0083] In the judgment rule, the ratio is set to 0.3. When the predicted horizontal offset is greater than the burial depth multiplied by 0.3, it indicates that the horizontal distance between the center of the radar anomaly and the center of the surface crack has exceeded one-third of the burial depth. At this time, it can be determined that the existing offset imaging has incorrectly repositioned the reflected energy of the tilted defect directly below the crack, that is, there is a radar offset repositioning error, and the actual defect and the surface crack are non-vertically related. Conversely, if the predicted horizontal offset is less than or equal to 0.3 times the burial depth, it is considered that the defect is basically located directly below the crack, and it is judged as vertically related.

[0084] To demonstrate the beneficial effects achieved in this step, a complete process from data acquisition to association determination is presented based on actual road detection data. The method is compared with the traditional vertical nearest association criterion, verifying its significant effectiveness in avoiding misjudgments.

[0085] A 100-meter-long section of a highway with karst development was inspected. The inspection vehicle traveled at 30 km / h, and for every pulse output by the ranging wheel (every 0.1 meters), a frame of surface grayscale image, a set of 16 elevation data points, and 8 radar waveforms were simultaneously acquired. At a distance of 54.2 meters, the surface image was filtered using anisotropic diffusion filtering and second-order partial derivative matrix eigenvalue ratio selection, revealing a transverse crack with a center transverse coordinate of 0.32 meters (0 point with the left edge of the road cross-section as the reference point), a crack width of 2.3 mm, and a crack orientation of 88 degrees (approximately perpendicular to the direction of travel). At the same distance, the internal radar data, after reverse time-shift imaging, identified a candidate anomaly on the depth profile with a center transverse coordinate of 0.31 meters, a burial depth of 1.2 meters, and a reflection intensity of 0.78 (normalized value).

[0086] According to the traditional vertical nearest correlation criterion, since the horizontal distance between the crack center (0.32 meters) and the anomaly center (0.31 meters) is only 0.01 meters, which is much less than 0.1 times the burial depth of 1.2 meters, it is determined to be vertically correlated. It is believed that the crack was caused by the cavity directly below, and it is recommended to drill and grout immediately.

[0087] The method described in this application is used for correlation judgment. First, seven measured input parameters are used: the transverse coordinate of the radar anomaly center is 0.31 meters, the burial depth is 1.2 meters, the reflection intensity is 0.78, the transverse coordinate of the surface crack center is 0.32 meters, the crack width is 2.3 millimeters, and the crack orientation is 88 degrees; these are used to form a measured input vector. Using a pre-constructed radial basis function interpolation mapping (this mapping is trained based on 208 sets of electromagnetic-mechanical coupled forward simulation datasets, with the width parameter taken as 0.5 times the maximum Euclidean distance of all sample points in the simulation dataset, and two sets of weighting coefficients obtained by solving a system of linear equations), the Euclidean distance between the measured vector and the 208 sample points is calculated, and then the radial basis function value is obtained, and the predicted value is obtained by weighted summation. The calculation output results are: the predicted defect tilt angle is 38 degrees, and the predicted horizontal offset is 0.78 meters. The burial depth is 1.2 meters, the set ratio is 0.3, and 0.3 times the burial depth is 0.36 meters. The predicted horizontal offset of 0.78 meters is greater than 0.36 meters, so it is determined to be a non-vertical correlation caused by radar offset repositioning error. That is, the actual defect is not located directly below the crack, but develops obliquely at an angle of 38 degrees, with a horizontal offset of 0.78 meters from the crack.

[0088] To verify the accuracy of the method, three-dimensional ground-penetrating radar cross-hole imaging was used to verify the location.

[0089] Two vertical holes, each 0.05 meters in diameter and 2.0 meters deep, were drilled at 0.3 meters and 1.1 meters on either side of the crack. The radar transmitting antenna was placed in the 0.3-meter hole, and the receiving antenna was placed in the 1.1-meter hole. The electromagnetic wave travel time was measured point by point. The inversion results showed that an ellipsoidal cavity exists at a depth of 1.15 to 1.35 meters, with its major axis at an angle of approximately 36 degrees to the horizontal plane. The center of the cavity is located 1.05 meters laterally (with the left edge as the 0 point), which is a horizontal offset of 0.73 meters relative to the crack center (0.32 meters).

[0090] The results are consistent with the defect inclination angle of 38 degrees and the horizontal offset of 0.78 meters predicted by the method of this invention, with errors within 0.05 meters and 2 degrees, respectively. Meanwhile, the cross-hole imaging at the location directly below the crack (0.31 meters laterally and 1.2 meters deep), which was misjudged by the traditional method, showed dense clay without any voids.

[0091] If grouting were performed according to the traditional method, construction workers would drill a hole directly below the crack (0.31 meters laterally) and inject grout at a pressure of 0.5 MPa, using 2.5 cubic meters of cement grout. However, because the actual cavity was located 1.05 meters laterally, the grout could not reach it, and the cavity remained. Three months later, during the rainy season, the road surface collapsed, creating a sinkhole approximately 1.2 meters in diameter and 1.5 meters deep. This caused the rear wheel of a light truck to become stuck, resulting in a two-hour traffic disruption.

[0092] After determining the defect to be non-vertical using the method of this invention, the maintenance unit, based on the predicted defect inclination angle of 38 degrees and horizontal offset of 0.78 meters, drilled and grouted at a horizontal distance of 1.05 meters, with a grouting pressure of 0.3 MPa and an injection volume of 1.8 cubic meters of cement grout, completely filling the cavity. Subsequent quarterly radar re-surveys over the following year showed no abnormal reflections, and no crack expansion or settlement occurred on the road surface.

[0093] This embodiment demonstrates that by generating a simulation dataset through pre-electromagnetic-mechanical coupling forward modeling and establishing a radial basis function nonlinear mapping, the present invention can accurately predict the tilt angle and horizontal offset of tilted defects, effectively identify radar offset repositioning errors under the traditional vertical nearest correlation criterion, avoid erroneous grouting and road collapse accidents caused by misjudgment, and significantly improve the reliability of joint detection of internal and surface defects of roads and the correctness of maintenance decisions.

[0094] Based on the judgment results, the association type between internal road defects and surface cracks is generated, and three types of states—vertical association, non-vertical association, or uncertain—are marked in the road 3D visualization data.

[0095] For each detection section (corresponding to a mileage coordinate), the association type at that section is generated based on the judgment result. The judgment result is derived by comparing the predicted horizontal offset with the burial depth: if the predicted horizontal offset is greater than 0.3 times the burial depth, the association type is marked as non-vertical association; if the predicted horizontal offset is less than or equal to 0.3 times the burial depth, the association type is marked as vertical association; if a clear judgment cannot be made due to low confidence (e.g., the confidence of the radial basis function interpolation output is less than 0.7), the association type is marked as uncertain. The confidence is calculated as follows: take the maximum value between the measured input vector and the radial basis function values ​​of all sample points in the simulated dataset. If this maximum value is greater than 0.7, the prediction is considered reliable; otherwise, it is considered uncertain. Each mileage coordinate and its corresponding association type, predicted defect inclination angle, predicted horizontal offset, radar anomaly burial depth, surface crack width, and other information are stored in a result data table, which is arranged in ascending order of mileage.

[0096] Construct 3D road visualization data. A 3D spatial grid is established, using the start and end mileage of the detected road segment as the longitudinal range, the road cross-section width (e.g., 7 meters) as the lateral range, and the detection depth (e.g., 3 meters) as the vertical range. The grid spacing is 0.1 meters longitudinally (consistent with the data acquisition cross-section interval), 0.05 meters laterally, and 0.05 meters vertically. Each grid point is assigned an initial color value: grid points above the road surface are set to gray, and grid points below the road surface without anomalies are set to semi-transparent blue.

[0097] The association type of each detection section is labeled. For sections marked as vertically associated, a cylindrical marker with a diameter of 0.2 meters is generated in the 3D visualization data, using the line connecting the coordinates of the surface crack center and the radar anomaly center as the axis. The cylinder extends from the road surface to the anomaly's burial depth, and its color is set to green. Simultaneously, a green dot with a radius of 0.05 meters is superimposed at the location of the road surface crack in this section. For sections marked as non-vertically associated, two independent markers are generated: one is a red dot at the surface crack center coordinates, and the other is a red sphere with a radius of 0.1 meters at the radar anomaly center coordinates (depth equal to burial depth). A red dashed line with a width of 0.02 meters connects the surface crack center and the anomaly center.

[0098] For sections marked as uncertain, a yellow dot is set at the center coordinates of the surface crack and a yellow sphere is set at the center coordinates of the radar anomaly, but no line is added between the two, and the text label "requires review" is added.

[0099] For road sections that are not vertically associated, the horizontal offset and inclination relationship between surface cracks and deep anomalies can be directly observed, which makes it easier for maintenance personnel to determine the precise grouting location.

[0100] Taking a specific road section as an example: at kilometer 54.2, it was determined to be non-vertically correlated. In the 3D view, a red dot appeared 0.32 meters laterally on the road surface, and a red sphere appeared at a depth of 1.2 meters and a width of 1.05 meters. Between the two, there was a red dashed line and the text label "Inclination angle 38 degrees, offset 0.78 meters". Based on this marking, maintenance personnel directly drilled and grouted at the width of 1.05 meters, avoiding incorrect construction. For uncertain cross-sections, cross-hole radar was used for verification before making a decision.

[0101] The above-mentioned visual annotations make the test results clear at a glance and effectively guide maintenance operations.

[0102] Example 2, please refer to Figure 2 As shown in this embodiment, a combined system for integrating ground-penetrating radar and road interior and surface defects includes: Data acquisition module: Simultaneously acquires road surface defect data and internal radar data along the road detection direction; Spatial registration module: Assigns unified mileage coordinates to each frame of surface defect data and each internal radar data, so that the same mileage position corresponds to one surface defect data unit and one internal radar data unit at the same time. Feature extraction module: performs image segmentation and edge detection on the surface defect data unit to extract the center coordinates, width and direction of surface cracks; performs offset imaging and anomaly recognition on the internal radar data unit to extract the center coordinates, burial depth and reflection intensity of radar anomalies; Forward modeling and mapping module: In advance, through the coupled forward modeling of electromagnetic wave propagation and solid mechanical response, a simulation dataset under multiple combinations of defect parameters is generated, and based on the simulation dataset, a nonlinear mapping function is established from the center coordinates, burial depth, reflection intensity of the radar anomaly, and the center coordinates, width, and orientation of the surface crack to the defect inclination angle and horizontal offset. The correlation judgment module inputs the center coordinates, width, and direction of the extracted surface cracks, as well as the center coordinates, burial depth, and reflection intensity of the radar anomaly into the nonlinear mapping function, and outputs the predicted defect tilt angle and the predicted horizontal offset. If the predicted horizontal offset is greater than the set ratio of the burial depth, it is determined to be a non-vertical correlation caused by the radar offset positioning error; otherwise, it is determined to be a vertical correlation. Visualization output module: Based on the judgment results, it generates the association type between internal road defects and surface cracks, and marks three types of states in the road 3D visualization data: vertical association, non-vertical association, or uncertain.

[0103] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for combining ground-penetrating radar with analysis of road interior and surface defects, characterized in that, include: Road surface defect data and internal radar data are collected synchronously along the road detection direction, and a unified mileage coordinate is assigned to each frame of road surface defect data and each internal radar data, so that the same mileage position corresponds to one surface defect data unit and one internal radar data unit at the same time. Image segmentation and edge detection are performed on the surface defect data unit to extract the center coordinates, width, and direction of the measured surface cracks; offset imaging and anomaly identification are performed on the internal radar data unit to extract the center coordinates, burial depth, and reflection intensity of the measured radar anomaly. A three-dimensional mesh array containing dielectric constant and elastic modulus values ​​is pre-constructed. The defect tilt angle and defect center coordinates of the defect region are set within this mesh array. On the same three-dimensional mesh array, the electromagnetic wave propagation equation is solved using the finite-difference time-domain method to obtain a simulated radar waveform sequence. The center coordinates, burial depth, and reflection intensity of the simulated radar anomaly are extracted from this sequence. A uniformly distributed standard reference load of 0.7 MPa is applied to the upper surface of the three-dimensional mesh array. The stress balance equation is solved using the finite element method to obtain the surface displacement field. The location with the largest vertical displacement gradient change in the surface displacement field is determined as the equivalent surface crack center coordinates. The width and orientation of the equivalent surface crack are determined based on the surface displacement field. The horizontal distance between the defect center coordinates and the equivalent surface crack center coordinates is calculated as the horizontal offset of this set of simulation data. The defect tilt angle, the horizontal offset, the center coordinates, burial depth and reflection intensity of the simulated radar anomaly, and the center coordinates, width and orientation of the equivalent surface crack are stored as simulated data records according to their corresponding relationships. The defect tilt angle and defect center coordinates are changed, and electromagnetic wave propagation calculation and solid mechanical response calculation are repeated to generate a simulated dataset composed of multiple sets of simulated data records. Based on the simulation dataset, a nonlinear mapping function is established using Gaussian radial basis functions, from the center coordinates, burial depth, and reflection intensity of the radar anomaly, as well as the center coordinates, width, and orientation of surface cracks, to the defect inclination angle and horizontal offset. The center coordinates, width, and orientation of the measured surface crack, as well as the center coordinates, burial depth, and reflection intensity of the measured radar anomaly, are input into the nonlinear mapping function. The function outputs the predicted defect inclination angle and the predicted horizontal offset. When the predicted horizontal offset is greater than 0.3 times the burial depth of the measured radar anomaly, the association type between the internal road defect and the measured surface crack is marked as non-vertical association. When the predicted horizontal offset is less than or equal to 0.3 times the burial depth of the measured radar anomaly, the association type is marked as vertical association. When the confidence level of the Gaussian radial basis function interpolation output is less than 0.7, the association type is marked as uncertain. Based on the association type, vertical association, non-vertical association, or uncertain state are marked in the road 3D visualization data.

2. The method for combining ground-penetrating radar with analysis of road interior and surface defects according to claim 1, characterized in that, The road surface defect data includes image information and elevation information of the road surface, and the internal radar data is electromagnetic wave reflection signal collected by a linear array composed of 8 transceiver radar antenna units; the surface defect data acquisition component and the internal radar data acquisition component are installed on the same detection vehicle and are synchronously triggered by the pulse signal of the same rotating coded ranging wheel.

3. The method for combining ground-penetrating radar with analysis of road interior and surface defects according to claim 1, characterized in that, Extracting the center coordinates, width, and orientation of the measured surface cracks, including: Perform anisotropic diffusion filtering on the grayscale image; The second-order partial derivative matrix at each pixel is calculated for the filtered grayscale image. Linear structure pixels are selected to form crack candidate regions based on the ratio of the two eigenvalues ​​of the second-order partial derivative matrix. Morphological skeleton extraction is performed on the candidate crack region to obtain a crack centerline with a single pixel width. The gray-level gradient is calculated along the normal direction of the crack centerline, and the gradient peak spacing is used as the crack width. The centerline of the crack is segmented and fitted with a straight line, and the direction angle of the fitted straight line is used as the crack direction.

4. The method for combining ground-penetrating radar with analysis of road interior and surface defects according to claim 1, characterized in that, Extracting the center coordinates, burial depth, and reflection intensity of the measured radar anomaly includes: Time zero-point correction and exponential gain compensation are performed on the eight waveform sequences in each internal radar data unit. The energy of each reflected wave is reversed along the propagation path and superimposed to generate a two-dimensional depth profile with depth as the vertical axis and lateral distance as the horizontal axis. The ratio of local energy to background noise energy is calculated point by point on the two-dimensional depth profile, and continuous pixel regions with a ratio greater than 4 are marked as candidate anomalous bodies. The lateral coordinates and depth values ​​of the geometric center of the candidate anomaly on the two-dimensional depth profile are calculated as the center coordinates and burial depth of the measured radar anomaly, and the maximum amplitude of all pixels in the candidate anomaly is taken as the reflection intensity.

5. The method for combining ground-penetrating radar with methods for detecting internal and surface defects of roads, as described in claim 1, is characterized in that... The three-dimensional mesh array has a horizontal length of 0.7 meters, a length in the road's forward direction of 0.1 meters, and a depth length of 3.0 meters, with a mesh spacing of 0.01 meters in all three directions. The dielectric constants of the surface layer, base layer, and subgrade are 6, 8, and 12, respectively, and the elastic moduli are 3000 MPa, 800 MPa, and 50 MPa, respectively. The defect region is an ellipsoid with a dielectric constant of 4 and an elastic modulus of 0.1 MPa. The defect tilt angle varies from 0 degrees to 60 degrees, and the defect center coordinates vary within the lateral range of the three-dimensional mesh array.

6. The method for combining ground-penetrating radar with analysis of road interior and surface defects according to claim 1, characterized in that, The nonlinear mapping function is established using the Gaussian radial basis function, including: The input and output parameters of each set of simulation data are extracted from the simulation dataset. The input parameters include the center coordinates, burial depth, reflection intensity of the simulated radar anomaly, and the center coordinates, width, and orientation of the equivalent surface crack. The output parameters include the defect tilt angle and horizontal offset. Using the input parameter vector of each sample point as the center, calculate the Euclidean distance between the input parameter vectors of any two sample points, and construct a radial basis function matrix of order number of sample points; Set 0.5 times the maximum value of the Euclidean distance between all sample points as the width parameter of the Gaussian radial basis function, solve the linear equation system formed by the radial basis function matrix and the output parameter matrix, and obtain the radial basis function weight coefficients corresponding to each sample point; For the measured input parameters, calculate the Euclidean distance between them and the input parameter vectors of each sample point, and perform a weighted summation based on the radial basis function weighting coefficients to obtain the predicted defect tilt angle and the predicted horizontal offset.

7. The method for combining ground-penetrating radar with analysis of road interior and surface defects according to claim 1, characterized in that, The confidence level is the maximum value of the Gaussian radial basis function values ​​corresponding to the measured input parameter vector and the input parameter vectors of each sample point in the simulated dataset.

8. The method for combining ground-penetrating radar with analysis of road interior and surface defects according to claim 1, characterized in that, When constructing the three-dimensional visualization data of the road, the starting and ending mileage of the detected road segment is used as the longitudinal range, the width of the road cross section is used as the lateral range, and the detection depth is used as the vertical range to establish a three-dimensional spatial grid. For vertical correlation, a green cylindrical marker is generated between the center of the measured surface crack and the center of the measured radar anomaly; for non-vertical correlation, a red dot is set at the center of the measured surface crack, and a red sphere is set at the predicted defect location and connected by a red dashed line; for uncertain states, yellow markers are set at the center of the measured surface crack and the center of the measured radar anomaly, and a verification label is added.

9. A combined system integrating ground-penetrating radar and road interior / surface defects, characterized in that, It includes a data acquisition module, a spatial registration module, a feature extraction module, a forward modeling and mapping construction module, an association judgment module, and a visualization output module; The data acquisition module is used to simultaneously acquire road surface defect data and internal radar data along the road detection direction; The spatial registration module is used to assign unified mileage coordinates to each frame of road surface defect data and each internal radar data. The feature extraction module is used to extract the center coordinates, width and orientation of the measured surface cracks, as well as the center coordinates, burial depth and reflection intensity of the measured radar anomaly. The forward modeling and mapping construction module is used to set the defect tilt angle and defect center coordinates of the defect region in a three-dimensional mesh array according to the method of fusion ground penetrating radar for road interior and surface defects as described in claim 1. Electromagnetic wave propagation calculation and solid mechanical response calculation under a uniformly distributed standard reference load of 0.7 MPa are performed respectively. The equivalent surface crack center coordinates are determined according to the surface displacement field, and the horizontal distance between the defect center coordinates and the equivalent surface crack center coordinates is calculated as the horizontal offset. A simulation dataset is generated and a Gaussian radial basis function nonlinear mapping function is established. The correlation judgment module is used to input the measured features into the nonlinear mapping function, and determine the vertical correlation, non-vertical correlation or uncertain state based on the comparison result between the predicted horizontal offset and 0.3 times the burial depth of the measured radar anomaly, as well as the confidence level. The visualization output module is used to annotate the vertical correlation, non-vertical correlation, or uncertain state in the road 3D visualization data.

Images