An underground pipeline surveying and mapping identification method combining ground penetrating radar and three-dimensional point cloud

Abstract

The application relates to the technical field of radar detection, and discloses a ground penetrating radar and three-dimensional point cloud combined underground pipeline surveying and mapping identification method, which comprises the following steps: acquiring ground penetrating radar, three-dimensional point cloud, attitude and mileage data, determining an underground reflection candidate area, calculating a local ground surface normal and arranging a local closed measurement unit; implementing bidirectional, double-normal-distance four-state observation, correcting and orthogonally decomposing reflection arrival time, obtaining four types of reflection components, completing trajectory synchronization and antenna coupling correction; establishing a common reflector observation group and a differential closed relationship based on updated common reflection components, determining a pipeline center position, a strike and a local equivalent radar wave propagation velocity, confirming or splitting the observation group according to basic residuals and parameter uncertainty, and generating a supplementary surveying action and output surveying and mapping results. Through separation of trajectory synchronization error and antenna coupling error and combination of differential closed inversion, the application improves the accuracy of underground pipeline positioning, strike identification and buried depth measurement.

Description

Technical Field

[0001] This invention relates to the field of radar detection technology, specifically to a method for mapping and identifying underground pipelines using a combination of ground-penetrating radar and three-dimensional point cloud. Background Technology

[0002] Currently, urban roads and work sites contain underground pipelines for water supply, drainage, gas, electricity, and communications. Some pipelines lack complete as-built documentation, and existing records may not accurately reflect their actual locations. Ground-penetrating radar (GPR) can identify pipelines based on reflected signals generated by changes in the electromagnetic parameters of the underground medium, while 3D point clouds can record surface morphology, feature boundaries, and the spatial location of measuring equipment. Combining underground reflection information with surface spatial information provides a data foundation for determining the center location, burial depth, and orientation of pipelines.

[0003] To address the aforementioned underground pipeline mapping needs, current operations typically involve controlling a ground-penetrating radar (GPR) to move along a predetermined survey line, continuously emitting electromagnetic pulses and receiving reflected signals from underground targets. Multiple radar data channels are then combined to form a radar profile. Processing personnel determine the target location based on the reflected waveform or hyperbolic characteristics and calculate the burial depth according to a preset radar wave propagation speed. The survey line trajectory is generally recorded by an odometer, satellite positioning equipment, or a total station. When using 3D point clouds, a surface model is primarily established by scanning the road surface, manhole covers, and surrounding features. The radar identification results are then converted to an engineering coordinate system to generate a plan view or 3D rendering of the underground pipelines.

[0004] Existing methods are heavily reliant on preset propagation velocities, and changes in the moisture content of the medium or backfill material will shift the calculated burial depth. Antenna distance from the ground, equipment attitude, and surface undulations also alter the reflection arrival time; these effects are often intertwined with underground propagation time, making them difficult to correct separately. Single-direction measurements provide limited constraints, making it difficult to detect trajectory misalignment, time shifts, and direction-related errors. Even with the addition of reverse survey lines, simply superimposing the results by position cannot utilize reciprocity relationships to separate common reflection components from directional disturbances. Point cloud data is primarily used for surface modeling and result registration, offering insufficient constraints on antenna phase center, local surface normal, and radar propagation path. Furthermore, pipeline center location, orientation, and local propagation velocity are parametrically coupled, resulting in multiple approximate solutions based solely on reflection data obtained from a single measurement height. Existing results are typically represented by identified points or fitted pipelines, lacking quantitative data such as closure residuals and parameter uncertainties. Even after local anomalies occur, it is difficult to accurately determine the remeasurement location, measurement direction, and acquisition status.

[0005] Therefore, this invention provides a method for mapping and identifying underground pipelines by combining ground penetrating radar and three-dimensional point cloud to address the shortcomings of existing technologies. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a ground-penetrating radar combined with 3D point cloud for underground pipeline mapping and identification. This method solves the problems in existing underground pipeline mapping processes caused by surface undulations, asynchronous equipment trajectories, variations in antenna normal distance, changes in antenna coupling state, uncertain propagation speed of underground media, and mixed reflections from multiple underground targets. These problems result in unstable reflection arrival time, large errors in determining pipeline location and burial depth, unreliable target association, and a lack of targeted supplementary surveys.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for mapping and identifying underground pipelines using ground-penetrating radar combined with 3D point cloud data, comprising the following steps: Acquire ground-penetrating radar data, 3D point cloud data, attitude data, and odometer data, and determine candidate areas for underground reflection; Calculate the local surface normal based on the 3D point cloud data within the candidate area and deploy local closed measurement units; Four-state observations of bidirectional and binormal distances were implemented in the measurement unit. Correct the reflection arrival time and orthogonally decompose to obtain the common reflection component, the directional odd component, the normal distance perturbation component, and the directional and normal distance interaction component; Trajectory synchronization and antenna coupling correction are performed based on each component, and the results are updated accordingly. A common reflector observation group is established using the updated common reflection component as the stable reflection arrival time; Establish a differential closure relationship for the observation group, determine the pipeline center position, direction and local equivalent radar wave propagation velocity, and obtain candidate propagation velocity values, basic residuals and parameter uncertainties; The observation group is confirmed or split based on the candidate propagation velocity values, each component, and the fundamental residual; Based on each component and the basic residual, supplementary measurement actions are generated, and the survey results are output according to the parameter uncertainty.

[0008] By adopting the above technical solution, the 3D point cloud provides a surface spatial reference and local surface normal, while the four-state observation with dual-direction and dual-normal distance provides independent conditions for direction and normal distance variation. The original reflection arrival times of candidate reflection events under the four observation states are corrected for air propagation time, fixed delay within the antenna, and orthogonal decomposition, respectively forming the common reflection component, directional odd component, normal distance perturbation component, and interactive component of direction and normal distance for underground targets. This allows for the separation of trajectory synchronization error, antenna coupling error, and stable reflection time of underground targets. Furthermore, by jointly obtaining pipeline spatial parameters and local equivalent radar wave propagation velocity through differential closure relationships, the stability of underground pipeline positioning, depth measurement, and target identification is improved.

[0009] Preferably, the acquisition of ground-penetrating radar data, 3D point cloud data, attitude data, and odometer data, and the determination of underground reflection candidate areas, includes: performing time-zero coarse correction on the ground-penetrating radar data; matching reflection events in adjacent radar channels whose reflection time difference and principal phase difference are not greater than preset trajectory time thresholds and trajectory phase thresholds, respectively, to form an initial reflection trajectory; determining the equipment trajectory based on the attitude data and odometer data, and projecting the initial reflection trajectory onto the surface represented by the 3D point cloud data according to the initial time index and the equipment trajectory; merging the initial reflection trajectories in adjacent survey lines whose projection position distance and reflection time difference are not greater than preset survey line merging distance thresholds and survey line merging time thresholds, respectively, to form underground reflection candidate areas.

[0010] By adopting the above technical solution, the temporal continuity, phase consistency, and surface projection position continuity of radar reflection events are combined to limit the subsequent closed measurement range and reduce the impact of isolated reflections, discontinuous reflections, and erroneous correlations between different survey lines on pipeline identification.

[0011] Preferably, the step of calculating the local surface normal based on the 3D point cloud data within the candidate area and deploying locally closed measurement units includes: Surface point clouds are extracted from candidate underground reflection regions, the covariance matrix of neighborhood point clouds is calculated, and the unit eigenvector corresponding to the smallest eigenvalue of the covariance matrix is ​​determined as the local surface normal. The initial extension direction of the candidate underground reflection region is fitted, and the projection of the initial extension direction onto the local tangent plane perpendicular to the local surface normal is determined as the longitudinal unit tangent, and the transverse unit tangent perpendicular to the longitudinal unit tangent is determined. Set up a central measuring point, and set up positive transverse measuring points, negative transverse measuring points, positive longitudinal measuring points, and negative longitudinal measuring points along the transverse unit tangent and the longitudinal unit tangent respectively to form a local closed measuring unit.

[0012] By adopting the above technical solution, the location of the measuring points is determined based on the actual surface geometry and the candidate reflection extension direction. This allows the lateral measuring points to be used to characterize the lateral position changes of the pipeline, and the longitudinal measuring points to be used to characterize the pipeline direction changes. Furthermore, the adjustment direction of the normal distance is made consistent with the local surface normal, reducing the horizontal projection offset caused by slopes and uneven surfaces.

[0013] Preferably, the implementation of bidirectional and bidirectional normal distance four-state observation in the measurement unit includes: using two different preset normal distances at the same ground projection position of each measuring point in the locally closed measurement unit, performing forward first normal distance measurement, reverse first normal distance measurement, forward second normal distance measurement, and reverse second normal distance measurement in sequence; recording the actual normal distance, attitude data, and mileage data of each state to determine the corresponding equipment trajectory; and resampling the ground penetrating radar data according to the three-dimensional point cloud coordinates so that the four states correspond to the same ground projection position.

[0014] By adopting the above technical solution, forward and reverse states form directional reciprocal observations, the first normal distance state and the second normal distance state form normal distance differential observations, and spatial resampling eliminates the differences in the surface projection positions between the four observation states, providing reflection arrival times with spatial correspondences for subsequent orthogonal decomposition.

[0015] Preferably, the step of correcting the reflection arrival time and orthogonally decomposing it to obtain the common reflection component, the directional odd component, the normal distance perturbation component, and the directional and normal distance interaction component includes: acquiring the extrinsic parameters between the pre-calibrated ground-penetrating radar antenna and the attitude sensor; determining the phase center of the transmitting antenna, the phase center of the receiving antenna, and the corresponding ground projection point based on the attitude data, the extrinsic parameters, and the local surface normal; calculating the air propagation time based on the distance from each phase center to the corresponding ground projection point and the air wave speed; extracting the reflection arrival time of candidate reflection events under four states from the ground-penetrating radar data as the original reflection arrival time; subtracting the air propagation time and the pre-calibrated fixed delay inside the antenna from the original reflection arrival time to obtain the corrected reflection arrival time; and performing orthogonal decomposition on the four corrected reflection arrival times of the same location and the same reflection event to obtain the common reflection component, the directional odd component, the normal distance perturbation component, and the directional and normal distance interaction component.

[0016] For a separate transmit / receive ground-penetrating radar, the air propagation path is the sum of the distance from the phase center of the transmitting antenna to the surface projection point of the transmitting end and the distance from the phase center of the receiving antenna to the surface projection point of the receiving end, and the air propagation time is the ratio of the sum of the distances to the air wave speed. For a transceiver integrated ground-penetrating radar, the air propagation path is twice the distance from the antenna phase center to the surface projection point, and the air propagation time is the ratio of the air propagation path to the air wave speed.

[0017] By employing the above technical solution, the reflection arrival times of candidate reflection events under four states are extracted from ground-penetrating radar data and determined as the original reflection arrival times. This clarifies the source of the reflection arrival times and their correspondence with the four-state observations. By subtracting the air propagation time and the internal fixed delay of the antenna and acquisition circuit from the original reflection arrival times, the air propagation time formed when the antenna phase center is above the ground surface is separated from the underground propagation time. This provides a unified underground propagation time reference for reflection arrival times under different ground slopes, different equipment attitudes, and different normal distances.

[0018] Specifically, the arrival times of four corrected reflections at the same location and for the same reflection event are orthogonally decomposed into directional and normal distance states. The common reflection component is the arithmetic mean of the arrival times of the four corrected reflections; the directional odd component is determined by averaging the time difference between forward and reverse measurements under two normal distance states; the normal distance perturbation component is determined by averaging the time difference between the second and first normal distance states in both forward and reverse measurements; and the directional and normal distance interaction component is determined by the change in the forward and reverse time difference between the two normal distance states.

[0019] By adopting the above technical solution, the common reflection component retains the underground reflection time that exists in all four observation states, the directional odd component represents the difference between forward and reverse trajectories and time synchronization, the normal distance disturbance component represents the change in normal distance and the change in antenna coupling, and the directional and normal distance interaction component represents the non-independent response between the change in direction and the change in normal distance, so that different error sources enter the corresponding correction, correlation and supplementary measurement processes respectively.

[0020] Preferably, the step of implementing trajectory synchronization and antenna coupling correction based on each component and updating it includes: constructing a synchronization correction model with the goal of minimizing the root mean square value of the directional odd component and the mileage ratio, mileage translation, and time synchronization offset of the reverse measurement as the optimization quantities; solving the synchronization correction model, redetermining the spatial position and time correspondence of the reverse measurement according to the obtained mileage ratio, mileage translation, and time synchronization offset, and re-performing spatial resampling and orthogonal decomposition.

[0021] By adopting the above technical solution, the directional odd component directly participates in the adjustment of the reverse measurement trajectory and time synchronization relationship, so that the forward and reverse data correspond to each other in spatial position and sampling time, reducing the reflection time difference caused by the slippage of the ranging wheel, the change of movement speed and the acquisition trigger offset.

[0022] Further preferably, the direct wave or stable surface reference reflection is extracted from a preset time window starting from time zero. Based on the normal distance disturbance component of the direct wave or stable surface reference reflection under different normal distance states, the system group delay corresponding to the change in normal distance is determined. The system group delay is deducted from the arrival time of each corrected reflection, and the common reflection component, directional odd component, normal distance disturbance component, and directional and normal distance interaction component are recalculated and updated.

[0023] By adopting the above technical solution, the system group delay caused by antenna coupling is calibrated using an early reference signal that does not contain changes in the propagation of underground targets or has small changes in the propagation of underground targets, thereby reducing the impact of changes in antenna off-ground state, ground contact state, and antenna operating state on underground reflection time.

[0024] Preferably, establishing a common reflector observation group based on the updated common reflection component as the stable reflection arrival time includes: retrieving reflection events where the stable reflection arrival time difference between adjacent measuring points is not greater than a preset continuity threshold; calculating the principal phase difference between reflection events, extracting the effective bandwidth from the spectrum and calculating the effective bandwidth relative difference; establishing association relationships for reflection events where the principal phase difference, effective bandwidth relative difference, and direction-normal distance interaction component are not higher than the associated phase difference threshold, bandwidth difference threshold, and direction-normal distance interaction threshold, respectively; and forming a common reflector observation group from the reflection events of at least three non-collinear measuring points.

[0025] The relative difference in effective bandwidth is the ratio of the absolute value of the difference in effective bandwidth between two reflection events to the average value of the two effective bandwidths. The associated phase difference threshold ranges from 15° to 60°, the bandwidth difference threshold ranges from 0.10 to 0.50, and the interaction threshold between direction and normal distance ranges from 0.03 ns to 1.00 ns.

[0026] By adopting the above technical solution, the stability of reflection time continuity, principal phase consistency, effective bandwidth consistency, and the interaction state of direction and normal distance are used together for reflection event correlation. It also requires that the observation group contain at least three non-collinear measurement points to reduce the probability that adjacent pipelines, intersecting pipelines, and blocky abnormal reflections are incorrectly classified into the same observation group.

[0027] Preferably, the step of establishing a differential closure relationship for the observation group, determining the pipeline center position, orientation, and local equivalent radar wave propagation velocity, and obtaining candidate propagation velocity values, fundamental residuals, and parameter uncertainties includes: representing the candidate pipeline segment axis using the center control point of the pipeline segment to be estimated and the orientation unit vector. The center control point is the point on the axis with the smallest distance between the candidate pipeline segment axis and the center of the local closure measurement unit, and the orientation unit vector is used to define the spatial extension direction of the candidate pipeline segment.

[0028] The underground propagation geometric distance is calculated based on the three-dimensional coordinates of each measuring point. For a separate transmit / receive ground penetrating radar, the underground propagation geometric distance is the sum of the shortest distance from the surface projection point of the transmitter to the axis of the candidate pipe segment and the shortest distance from the surface projection point of the receiver to the axis of the candidate pipe segment; for an integrated transmit / receive ground penetrating radar, the underground propagation geometric distance is twice the shortest distance from the surface projection point of the antenna to the axis of the candidate pipe segment.

[0029] A differential propagation time relationship is established based on the difference in stable reflection arrival times between different measuring points. The differential propagation time is determined jointly by the difference in underground propagation geometric distances between corresponding measuring points and the local equivalent radar wave propagation velocity. The common system time delay and common pipeline radius response within the same common reflector observation group are mutually canceled after subtracting the stable reflection arrival times.

[0030] A transverse antisymmetric time difference is generated based on the difference in stable reflection arrival times between the positive and negative transverse measuring points. A longitudinal antisymmetric time difference is generated based on the difference in stable reflection arrival times between the positive and negative longitudinal measuring points. A transverse second-order closure variable is generated based on the difference between the sum of the stable reflection arrival times of the positive and negative transverse measuring points and twice the stable reflection arrival time of the central measuring point. A longitudinal second-order closure variable is generated based on the difference between the sum of the stable reflection arrival times of the positive and negative longitudinal measuring points and twice the stable reflection arrival time of the central measuring point. Candidate propagation velocity values ​​are obtained from these difference closure relationships.

[0031] By adopting the above technical solution, the differential propagation time is used to constrain the relationship between the propagation geometric distance and the local equivalent radar wave propagation speed, the lateral antisymmetric time difference is used to constrain the lateral position of the pipeline, the longitudinal antisymmetric time difference is used to constrain the pipeline direction, and the lateral and longitudinal second-order closure quantities are used to constrain the pipeline burial depth and propagation speed, so that the pipeline position, direction, burial depth and propagation speed can be jointly determined within the same local measurement unit.

[0032] Further preferredly, a joint objective function is constructed based on the differential propagation residual, antisymmetric closure residual, second-order closure residual, and the interaction components of direction and normal distance, and preset weights are assigned to each type of residual; a weighted robust optimization method is used to solve for the center control point, direction unit vector, and local equivalent radar wave propagation velocity of the pipeline segment to be estimated, thereby determining the pipeline center position, direction, and local equivalent radar wave propagation velocity; the differential propagation residual, antisymmetric closure residual, and second-order closure residual are determined as the basic residuals, and the parameter covariance matrix is ​​determined based on the sensitivity relationship of the residuals to the parameters to be estimated, the observation weights, and the residual variance, and the parameter uncertainty is determined by the parameter covariance matrix.

[0033] By adopting the above technical solution, the joint objective function simultaneously utilizes the propagation time difference, the closure relationship of symmetrical measuring points, and the interaction state between direction and normal distance, and reduces the impact of anomalous reflection and mixed reflection on parameter solving through robust weights. The parameter covariance matrix is ​​used to characterize the uncertainty of pipeline plane position, burial depth, orientation, and local equivalent radar wave propagation velocity.

[0034] Preferably, the step of confirming or splitting the observation group based on the candidate propagation velocity values, each component, and the fundamental residuals includes: obtaining candidate propagation velocity values ​​from different valid observation pairs; determining the ratio of the standard deviation to the mean of the candidate propagation velocity values ​​as the propagation velocity closure index; and determining the root mean square value of the interaction components of the direction and normal distances within the common reflector observation group as the interaction closure index. When the propagation velocity closure index, the interaction closure index, and each fundamental residual are all not higher than their corresponding thresholds, the common reflector observation group is confirmed as a single linear pipeline.

[0035] When any index or any basic residual exceeds the corresponding threshold, subgroups are formed according to the principal phase of the reflection event, the arrival time of stable reflection, the spatial projection position, and the interaction components of the direction and normal distance. Differential closure relationships are established for each subgroup. If no subgroup meets the aforementioned conditions after reaching the preset maximum number of retests, the absolute value of the second-order closure quantity obtained by combining reflection events from different subgroups is determined as the mixed second-order closure residual. When the mixed second-order closure residual exceeds the diagonal closure threshold, the corresponding target is marked as a nonlinear underground target group.

[0036] The propagation speed closure threshold ranges from 0.03 to 0.20, and the threshold for the interactive closure index is the interaction threshold between the direction and the normal distance. The diagonal closure threshold is not less than 0.05 ns.

[0037] By adopting the above technical solutions, the propagation velocity closure index evaluates whether the same observation group can be explained by a unified local propagation velocity, the interactive closure index evaluates whether the changes in direction and normal distance can be independently separated, and the basic residual evaluates the consistency between the linear propagation model and the measured data, thereby distinguishing between single pipelines, mixed reflections of multiple pipelines, and nonlinear underground targets.

[0038] Preferably, the step of generating supplementary measurement actions based on each component and the basic residual, and outputting the mapping results according to the parameter uncertainty, includes: determining the root mean square value of the directional odd component as the directional residual, determining the root mean square value of the normal distance disturbance component as the normal distance residual, determining the root mean square value of the directional and normal distance interaction component as the directional and normal distance interaction residual, determining the root mean square value of the differential propagation residual as the differential propagation residual index, and determining the larger of the root mean square value of the antisymmetric closed residual and the root mean square value of the second-order closed residual as the closed residual.

[0039] When the directional residual exceeds the corresponding threshold, a reverse remeasurement action on the same trajectory is generated; when the normal distance residual exceeds the corresponding threshold, a supplementary measurement action to increase the normal distance state is generated; when the interaction residual between directional and normal distance exceeds the corresponding threshold, a supplementary measurement action to add diagonal measuring points and split the common reflector observation group is generated; when the differential propagation residual index exceeds the corresponding threshold, a supplementary measurement action to adjust the center of the locally closed measurement unit is generated; when the closure residual exceeds the corresponding threshold, a supplementary measurement action to change the measurement baseline is generated. After the supplementary measurement is completed, four-state observations are re-implemented and differential closure relationships are established.

[0040] By adopting the above technical solutions, various residuals correspond to trajectory synchronization state, antenna coupling state, reflection correlation state, measurement unit position state, and propagation parameter identifiable state, respectively. This allows the supplementary measurement position, supplementary measurement direction, normal distance state, and measurement baseline to be determined based on the source of error, reducing data redundancy caused by indiscriminate supplementary measurements.

[0041] Further optimized, the standard deviation of pipeline plane position and the standard deviation of pipeline burial depth are used to characterize parameter uncertainty. When each closure index, each residual, the standard deviation of plane position, and the standard deviation of burial depth do not exceed the corresponding threshold, the axes of adjacent candidate pipe segments whose distance between adjacent endpoints is not greater than the preset connection distance threshold and whose directional angle is not greater than the preset connection angle threshold are connected to form the three-dimensional centerline of the underground pipeline; if the closure condition or accuracy condition is still not met after reaching the preset maximum number of retests, an unclosed anomaly mark is output.

[0042] By adopting the above technical solutions, the surveying results are simultaneously constrained by closure consistency, residual level and parameter uncertainty, avoiding the connection of underground pipelines based solely on reflection intensity or single location results, and preserving the abnormal state of underground targets that cannot be stably interpreted.

[0043] This invention provides a method for mapping and identifying underground pipelines using a combination of ground-penetrating radar and 3D point cloud analysis. It offers the following advantages: 1. This invention utilizes three-dimensional point clouds to determine the local surface normal within the candidate area of ​​underground reflection, and accordingly deploys local closed measurement units and controls the normal distance state of ground penetrating radar, so that radar data from different measuring points and under different normal distance states have a unified surface spatial reference, which can reduce the impact of surface slope, equipment attitude changes and antenna surface projection offset on the location results of underground pipelines.

[0044] 2. This invention uses four-state observations of two-way and two-way normal distances to orthogonally decompose the corrected reflection arrival time, obtaining the common reflection component, directional odd component, normal distance disturbance component, and directional and normal distance interaction component. The corresponding components are then used to implement trajectory synchronization correction and antenna coupling correction, thereby separating forward and reverse trajectory errors, time synchronization errors, and antenna coupling changes, and improving the accuracy of stable reflection arrival time.

[0045] 3. This invention establishes a differential closure relationship based on a common reflector observation group, jointly determines the pipeline center position, direction and local equivalent radar wave propagation velocity, and confirms the observation group, splits the target and performs directional supplementary measurement based on the propagation velocity closure index, basic residual and parameter uncertainty. It can reduce the impact of local medium propagation velocity changes and multi-target mixed reflection on the mapping results, and output anomaly marks for areas that do not meet the closure conditions or accuracy conditions. Attached Figure Description

[0046] Figure 1 This is a framework diagram of the underground pipeline mapping and identification system of the present invention; Figure 2 This is a flowchart illustrating the overall process of the underground pipeline mapping and identification method of the present invention. Figure 3 This is a flowchart of the joint preliminary test and candidate region determination process of the present invention; Figure 4 This is a flowchart of the point cloud modeling and local closure measurement unit layout of the present invention; Figure 5 This is a flowchart of the four-state reciprocal acquisition and error correction process of the present invention; Figure 6 This is a flowchart of the differential closure inversion and pipeline identification process of the present invention; Figure 7 This is a flowchart of the residual assessment and targeted supplementary testing process of the present invention; Figure 8 This is a graph showing the changes in four-state reflection time and inversion residuals according to the present invention. Figure 9 This is a comparison curve of the absolute error of burial depth for different underground pipeline mapping methods of the present invention. Detailed Implementation

[0047] The technical solutions in 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, and 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.

[0048] In this embodiment, the unit of ground-penetrating radar reflection arrival time is ns, the unit of spatial coordinates, measurement distance, pipeline burial depth, and measurement baseline is meters, and the unit of radar wave propagation speed is m / ns. Ground-penetrating radar data, 3D point cloud data, attitude data, mileage data, and normal displacement data are recorded using a unified time reference and converted to the same engineering coordinate system.

[0049] In this embodiment, the stable reflection arrival time refers to the reflection arrival time represented by the updated common reflection component after completing air propagation time correction, antenna internal fixed delay correction, trajectory synchronization correction, and antenna coupling correction.

[0050] In this embodiment, the basic residuals include differential propagation residuals, antisymmetric closure residuals, and second-order closure residuals. The antisymmetric closure residuals include transverse antisymmetric closure residuals and longitudinal antisymmetric closure residuals, while the second-order closure residuals include transverse second-order closure residuals and longitudinal second-order closure residuals.

[0051] See attached document Figure 1 This invention provides an underground pipeline mapping and identification system that combines ground-penetrating radar and three-dimensional point cloud, including a joint acquisition module, a point cloud modeling module, a reciprocal correction module, a closed inversion module, and a result evaluation module.

[0052] The joint acquisition module establishes data connections with the point cloud modeling module and the reciprocal calibration module respectively. The point cloud modeling module establishes data connections with the reciprocal calibration module and the closed inversion module respectively. The reciprocal calibration module establishes a data connection with the closed inversion module. The closed inversion module establishes a data connection with the result evaluation module. The result evaluation module establishes a supplementary measurement control connection with the joint acquisition module.

[0053] The joint acquisition module transmits 3D point cloud data, attitude data, equipment reference point trajectory, and sampling time labels to the point cloud modeling module, and ground penetrating radar data, odometer data, normal displacement status, and sampling time labels to the reciprocal calibration module. The point cloud modeling module transmits local surface normal, antenna phase center coordinates, antenna-to-surface distance, and forward and reverse measurement position correspondence to the reciprocal calibration module, and transmits local tangent planes, local closed measurement unit coordinates, surface boundaries, and feature boundaries to the closed inversion module.

[0054] The reciprocity correction module transmits the stable reflection arrival time after air propagation time correction, trajectory synchronization correction, and antenna coupling correction to the closed-loop inversion module. It also transmits the common reflection component, directional odd component, normal distance perturbation component, and the interaction component between directional and normal distance. The closed-loop inversion module transmits the underground pipeline center location, pipeline direction, local equivalent radar wave propagation velocity, candidate propagation velocity values, basic residuals, and parameter uncertainties to the results evaluation module. Based on each component, basic residual, closure index, and parameter uncertainty, the results evaluation module generates supplementary measurement instructions and transmits the supplementary measurement location, supplementary measurement direction, normal distance status, and measurement baseline to the joint acquisition module.

[0055] The joint acquisition module includes a ground-penetrating radar acquisition unit, a 3D point cloud acquisition unit, an attitude acquisition unit, a mileage acquisition unit, a time synchronization unit, a normal displacement unit, a rigid mounting frame, and a field control terminal.

[0056] The ground-penetrating radar (GPR) acquisition unit includes a transmitting antenna, a receiving antenna, a radar controller, and a data acquisition unit. The transmitting antenna transmits electromagnetic pulses into the underground medium, and the receiving antenna receives reflected signals generated by abrupt changes in the dielectric parameters of the underground environment. The GPR acquisition unit generates radar A-scan data according to a set spatial sampling interval or temporal sampling interval, and assembles radar B-scan data according to the sequence of equipment movement.

[0057] The ground-penetrating radar acquisition unit adopts either an integrated transceiver structure or a separate transceiver structure. When using a separate transceiver structure, the phase centers of the transmitting and receiving antennas are spatially calibrated separately. When using an integrated transceiver structure, the phase centers of the transmitting and receiving antennas share the same calibration location.

[0058] The 3D point cloud acquisition unit includes a 3D laser scanning device, a point cloud acquisition controller, and a point cloud storage device. The 3D laser scanning device is used to acquire spatial point clouds of road surfaces, hardened ground, soil surfaces, manhole covers, valve boxes, inspection wells, curbs, walls, guardrails, and other ground objects. The 3D point cloud acquisition unit records the sampling time of each point cloud frame and converts the point cloud frames to the coordinate system of the joint acquisition equipment.

[0059] The attitude acquisition unit is mounted on a rigid frame and is used to acquire the heading, pitch, and roll states of the combined acquisition equipment. The sampling frequency of the attitude acquisition unit is no less than the spatial position update frequency of the ground-penetrating radar (GPR), ensuring that each radar A-scan data has corresponding attitude data. The attitude data output by the attitude acquisition unit is used to transform the phase center of the GPR antenna from the equipment coordinate system to the engineering coordinate system.

[0060] The odometer acquisition unit is mounted on the traveling mechanism of the combined acquisition equipment and is used to record the initial travel distance of the equipment along the survey line. The output of the odometer acquisition unit serves as the initial sequential position of the radar A-scan data, not the final spatial position. The final equipment trajectory is determined jointly by 3D point cloud data, attitude data, odometer data, and the synchronous correction results of forward and reverse trajectories.

[0061] The time synchronization unit is connected to the ground-penetrating radar acquisition unit, the 3D point cloud acquisition unit, the attitude acquisition unit, the odometer acquisition unit, and the normal displacement unit. The time synchronization unit assigns a unified time reference to each acquisition unit and records the trigger delay of each unit. The initial time correspondence generated by the time synchronization unit is used to establish a data index, and the remaining time synchronization offset between forward and reverse measurements is further calculated by the reciprocity correction module.

[0062] The normal displacement unit includes a guiding mechanism, a driving mechanism, a displacement detection mechanism, and a locking mechanism. The direction of movement of the guiding mechanism is adjusted according to the local surface normal output by the point cloud modeling module. The driving mechanism controls the movement of the ground-penetrating radar antenna phase center between the first normal distance state and the second normal distance state. The normal distance is the distance from the ground-penetrating radar antenna phase center along the local surface normal to the corresponding surface projection point. The displacement detection mechanism records the actual normal distance of the ground-penetrating radar antenna phase center, and the locking mechanism keeps the position of the ground-penetrating radar antenna phase center unchanged during radar acquisition.

[0063] The distance difference between the first normal range state and the second normal range state is determined based on the center frequency of the ground-penetrating radar, antenna structure, ground coupling conditions, and target burial depth. Preferably, the distance difference is 0.02m to 0.20m. Both the first and second normal range states meet the effective working distance requirements of the ground-penetrating radar antenna, and the signal-to-noise ratio of the underground target reflection in the second normal range state is not lower than a preset signal-to-noise ratio threshold.

[0064] A rigid mounting frame is used to maintain a fixed positional relationship between the ground-penetrating radar acquisition unit, the 3D point cloud acquisition unit, and the attitude acquisition unit. The rigid mounting frame is equipped with at least three point cloud identification markers, and the spatial distance between each marker is pre-calibrated. Before and after each field measurement, the 3D point cloud acquisition unit acquires the spatial coordinates of the point cloud identification markers. If the difference between the measured distance and the calibrated distance between the point cloud identification markers exceeds the allowable deviation of the external parameters, the original external parameters are discontinued and recalibrated. The allowable deviation of the external parameters is preferably 0.003m to 0.020m.

[0065] The point cloud modeling module includes a point cloud preprocessing unit, a surface extraction unit, a surface normal calculation unit, a ground feature identification unit, and a closure unit planning unit.

[0066] The point cloud preprocessing unit receives the 3D point cloud data output by the joint acquisition module and performs time alignment, coordinate transformation, outlier removal, duplicate point merging, and dynamic object point removal. For forward and reverse point clouds of the same area, the point cloud preprocessing unit performs spatial registration by fixing the ground features and road surface geometry, so that the forward and reverse point clouds are in the same engineering coordinate system.

[0067] The surface extraction unit extracts road surfaces, hardened ground surfaces, and soil surfaces as surface point clouds based on point cloud elevation continuity, local slope, surface normal variation, and spatial connectivity. Point clouds at the top of manhole covers and valve boxes are retained as surface point clouds when their elevations are continuous with the surrounding road surface. The edges of manhole covers, valve boxes, manholes, curbs, and potholes are recorded as abrupt surface boundaries.

[0068] The surface normal calculation unit establishes a point cloud neighborhood around each radar sampling location or locally closed measurement location, and fits the local surface plane using the neighborhood point cloud. The selection of the point cloud neighborhood simultaneously satisfies the requirements for the number of effective points and the consistency of the local plane. Preferably, the horizontal radius of the point cloud neighborhood is 0.15m to 0.80m, and the number of effective points is 30 to 300.

[0069] The surface normal calculation unit determines the direction with the least dispersion of the neighboring point cloud as the local surface normal, and adjusts the direction of the local surface normal to above the surface. The surface normal calculation unit establishes two mutually orthogonal local tangents within a local tangent plane perpendicular to the local surface normal. One local tangent is set as the longitudinal unit tangent based on the initial extension direction of the underground reflection candidate region, and the other local tangent is set as the transverse unit tangent perpendicular to the longitudinal unit tangent.

[0070] The surface normal calculation unit determines the surface roughness based on the characteristic value distribution of the neighboring point cloud. When the surface roughness exceeds the roughness threshold, the neighborhood radius is adjusted and refitted. Locations that still exceed the roughness threshold after refitting are not used as locations for dual normal distance acquisition, and alternative locations are selected by the closed-cell planning unit from adjacent passable areas. The roughness threshold is determined based on the point cloud statistical results of the smooth road reference area, preferably the sum of the mean roughness of the reference area and 2 to 4 times the standard deviation.

[0071] The feature identification unit identifies manhole covers, valve boxes, inspection wells, curbs, walls, guardrails, vehicles, and green facilities based on point cloud geometry, elevation changes, and spatial connectivity. Dynamic point clouds corresponding to vehicles and pedestrians are not included in the surface fitting. The boundaries of manhole covers, valve boxes, and inspection wells serve as auxiliary locations for pipeline topological connections. Curbstones, walls, guardrails, vehicles, and green facilities serve as access restriction boundaries for locally closed measurement units.

[0072] The closed-loop planning unit establishes a local closed measurement unit based on the initial location, initial extension direction, local surface normal, and passable area of ​​the candidate underground reflection area. The local closed measurement unit includes a central measuring point, a positive transverse measuring point, a negative transverse measuring point, a positive longitudinal measuring point, and a negative longitudinal measuring point. When there are cross reflections, mixed reflections from multiple targets, or significant directional inconsistencies, the local closed measurement unit also includes four diagonal measuring points.

[0073] The central measuring point is located near the initial surface projection of the candidate underground reflection location. Positive and negative transverse measuring points are symmetrical about the central measuring point, as are positive and negative longitudinal measuring points. The transverse measurement baseline is preferably 0.2 to 1.2 times the expected pipeline burial depth, and the longitudinal measurement baseline is preferably 0.5 to 2 times the transverse measurement baseline.

[0074] The closure unit planning unit performs a accessibility check on the predetermined measurement points. When the predetermined measurement point is within the obstacle range, two symmetrical measurement points are moved synchronously along the corresponding local tangent to ensure that the adjusted two measurement points still maintain a symmetrical relationship with respect to the central measurement point. If the symmetrical relationship cannot be maintained, the closure unit planning unit records the actual measurement point coordinates, and the closure inversion module constructs the propagation geometry according to the actual measurement point coordinates, without using the nominal baseline to replace the actual baseline.

[0075] The reciprocity correction module includes a four-state control unit, an air propagation correction unit, a reciprocity decomposition unit, a trajectory synchronization correction unit, and an antenna coupling correction unit.

[0076] The four-state control unit receives the coordinates of the measuring points, the local surface normal, and the normal distance status from the local closed measuring unit. It then controls the joint acquisition module to sequentially perform forward first normal distance measurement, reverse first normal distance measurement, forward second normal distance measurement, and reverse second normal distance measurement at the same surface projection position of each measuring point.

[0077] The four-state control unit ensures that forward and reverse measurements cover the same set of surface spatial locations. After the forward and reverse measurements are completed, the reciprocity correction module resamples the ground-penetrating radar data based on the three-dimensional point cloud coordinates, so that the radar reflection data in the four states correspond to the same surface projection location.

[0078] The air propagation correction unit determines the phase center of the transmitting antenna, the phase center of the receiving antenna, the ground projection point of the transmitting end, and the ground projection point of the receiving end based on the attitude data, sensor extrinsic parameters, local ground normal, and actual normal distance under each observation condition. The air propagation correction unit subtracts the air propagation time between the antenna phase center and the ground projection point from the original reflection arrival time, and also subtracts the internal fixed delay of the antenna and acquisition circuit.

[0079] The reciprocal decomposition unit orthogonally decomposes the arrival times of four corrected reflections at the same measurement location and for the same candidate reflection event, obtaining common reflection components, directional odd components, normal range perturbation components, and directional-normal range interaction components. The common reflection component is used for subsequent underground propagation inversion. The directional odd components are used to evaluate the consistency of forward and reverse device trajectories and time synchronization. The normal range perturbation components are used to evaluate the antenna's off-ground state and antenna coupling state. The directional-normal range interaction components are used to evaluate whether there is a non-independent response between changes in motion direction and changes in normal range.

[0080] The trajectory synchronization correction unit constructs a synchronization correction model. The synchronization correction model aims to minimize the root mean square value of the directional odd components of valid candidate reflection events, with the mileage ratio of reverse measurements, mileage translation, and time synchronization offset as the optimization parameters. After solving the synchronization correction model, spatial resampling is performed on both forward and reverse measurements, and a four-state orthogonal decomposition is re-executed.

[0081] The antenna coupling correction unit extracts the direct wave or stable ground reference reflection from a preset time window starting at time zero. It uses the normal distance perturbation component of the reference reflection to determine the time offset corresponding to the change in normal distance and defines the time offset as the system group delay. The antenna coupling correction unit subtracts the system group delay from the arrival time of each corrected reflection and re-executes the four-state orthogonal decomposition.

[0082] The directional residual threshold, normal distance residual threshold, and directional-normal distance interaction threshold are determined by the field reference area. The reference area is selected from regions with stable surface structures, continuous point cloud normals, and no significant underground target overlap. Each threshold is calculated as the sum of the mean of the corresponding reference residual and 2 to 4 times the standard deviation. The preferred directional residual threshold is no less than 0.05 ns and no more than 1.50 ns, the preferred normal distance residual threshold is no less than 0.05 ns and no more than 2.00 ns, and the preferred directional-normal distance interaction threshold is no less than 0.03 ns and no more than 1.00 ns.

[0083] The closed inversion module includes reflection correlation unit, pipeline modeling unit, difference calculation unit, parameter solving unit, and target splitting unit.

[0084] The reflection correlation unit receives the common reflection component, the directional odd component, the normal distance disturbance component, the directional and normal distance interaction component, the reflection amplitude, the reflection principal phase, the effective reflection bandwidth, and the coordinates of the measurement point. Based on the time continuity, spatial continuity, phase consistency, bandwidth consistency, and interaction residuals, it establishes a common reflector observation group.

[0085] The same candidate reflection event can simultaneously enter multiple initial common reflector observation groups. After closure inversion is completed, the reflection association unit determines its final assignment based on the fundamental residuals, propagation velocity closure index, and interaction closure index of each observation group. Reflection events that cannot satisfy the closure conditions of any observation group are marked as unidentified reflection events.

[0086] The pipeline modeling unit represents each local candidate pipe segment as a spatial linear target. A local pipe segment includes a central control point and a direction unit vector. The central control point is set as the closest point from the pipeline axis to the center of the local closed measurement unit to eliminate parameter redundancy caused by arbitrary movement of the central control point along the pipeline axis.

[0087] The differential calculation unit calculates the underground propagation geometric distance from each measuring point to the candidate pipeline axis based on the actual three-dimensional coordinates of the measuring points, and calculates the steady-state reflection arrival time difference between different measuring points. The differential calculation unit also generates the lateral antisymmetric time difference, the longitudinal antisymmetric time difference, the lateral second-order closure quantity, and the longitudinal second-order closure quantity, and calculates the candidate propagation velocity values ​​based on different differential closure relationships.

[0088] The lateral antisymmetric time difference is used to constrain the projected position of the pipeline axis in the local lateral direction, while the longitudinal antisymmetric time difference is used to constrain the deviation between the pipeline orientation and the initial extension direction. The lateral second-order closure and the longitudinal second-order closure are used to constrain the local pipeline burial depth and the local equivalent radar wave propagation velocity.

[0089] The parameter solving unit employs a weighted robust optimization approach to jointly determine the pipeline center location, pipeline direction, and local equivalent radar wave propagation velocity. During the iteration process, the parameter solving unit updates the observation weights based on the residual magnitude. Observations with large directional residuals, normal distance residuals, directional and normal distance interaction residuals, or large fundamental residuals have their weights reduced, while observations with stable reflection amplitudes, consistent phases, and continuous bandwidth have their weights increased to minimize the impact of anomalous observations on the solution results.

[0090] The target splitting unit determines whether a common reflector observation group corresponds to a single pipeline based on the propagation velocity closure index, the interaction closure index, the fundamental residual, and the mixed second-order closure residual. When the common reflector observation group cannot be explained by a single linear target and a single local equivalent radar wave propagation velocity, the target splitting unit divides the common reflector observation group into multiple subgroups based on the spatial projection position, the principal phase of the reflection, the stable reflection arrival time, and the interaction components of direction and normal distance, and performs differential closure inversion on each subgroup separately.

[0091] The results evaluation module includes a residual classification unit, an accuracy evaluation unit, a supplementary measurement planning unit, a centerline connection unit, and a results output unit.

[0092] The residual classification unit calculates the directional residual, normal distance residual, directional and normal distance interactive residual, differential propagation residual index, antisymmetric closure residual, second-order closure residual, and closure residual, respectively. Each type of residual corresponds to the trajectory synchronization state, antenna coupling state, reflection correlation state, measurement unit center position state, burial depth propagation velocity identifiable state, and propagation model consistency state, respectively.

[0093] The accuracy assessment unit calculates the parameter covariance matrix based on the sensitivity relationship of the parameter to be determined, the observation weights, and the residual variance. The accuracy assessment results include the standard deviations of pipeline horizontal position, pipeline elevation, pipeline burial depth, pipeline alignment, and the local equivalent radar wave propagation velocity. At least the standard deviations of horizontal position and burial depth are used to characterize the parameter uncertainty.

[0094] The supplementary measurement planning unit generates restricted candidate supplementary measurement actions based on each component and the basic residual. When the direction residual exceeds the direction residual threshold, candidate supplementary measurement actions include reverse remeasurement of the same trajectory and time-synchronized verification. When the normal distance residual exceeds the normal distance residual threshold, candidate supplementary measurement actions include adding a normal distance state and repeating the baseline normal distance measurement. When the interaction residual between direction and normal distance exceeds the interaction threshold between direction and normal distance, candidate supplementary measurement actions include adding diagonal measurement points and splitting the common reflector observation group.

[0095] When the differential propagation residual index exceeds the differential propagation residual threshold, the supplementary measurement planning unit adjusts the center of the local closed measurement unit according to the sign and spatial distribution of the differential propagation residual in the lateral and longitudinal directions. When the closed residual exceeds the closed residual threshold, the supplementary measurement planning unit changes the lateral or longitudinal measurement baseline, or adds a measurement line with a larger angle to the current pipeline direction.

[0096] Centerline connection units connect the three-dimensional centerlines of underground pipelines based on the distance between adjacent endpoints of adjacent candidate pipe segments, the included angle of orientation, burial depth, propagation velocity, and ground feature boundaries. Manholes, valve boxes, and well chamber boundaries serve as candidate topological boundaries and are not directly used as the center locations of underground pipelines. A corresponding topological relationship is established only when radar observation continuity, spatial continuity, and well chamber boundaries simultaneously meet the connection conditions.

[0097] The output unit outputs the three-dimensional centerline of the underground pipeline, the pipeline's plane coordinates, the pipeline's center elevation, the pipeline's burial depth, the local equivalent radar wave propagation velocity, the spatial confidence range, the quality grade, the original radar observation index, the supplementary measurement records, and the unclosed anomaly markers.

[0098] See attached document Figure 2 This invention provides a method for mapping and identifying underground pipelines using a combination of ground-penetrating radar and 3D point cloud. The method operates based on the aforementioned underground pipeline mapping and identification system and includes steps S100, S200, S300, S400, and S500.

[0099] Step S100 involves acquiring ground-penetrating radar data, 3D point cloud data, attitude data, and odometer data, and identifying candidate areas for underground reflection. This is performed through the joint acquisition module.

[0100] Step S200: Calculate the local surface normal based on the 3D point cloud data within the candidate area of ​​underground reflection, establish the local surface coordinates, and deploy local closed measurement units. This is executed through the point cloud modeling module.

[0101] Step S300 involves performing bidirectional and bidirectional normal range four-state observations within a locally closed measurement unit, correcting the reflection arrival time and performing orthogonal decomposition, completing trajectory synchronization correction and antenna coupling correction, and updating the common reflection component, directional odd component, normal range perturbation component, and directional and normal range interaction component. This is executed through the reciprocal correction module.

[0102] Step S400: Using the updated common reflection component as the stable reflection arrival time, establish a common reflector observation group and differential closure relationship, determine the center location, direction, and local equivalent radar wave propagation velocity of the underground pipeline, obtain candidate propagation velocity values, basic residuals, and parameter uncertainties, and confirm or split the common reflector observation group. This is executed through the closure inversion module.

[0103] Step S500: Determine the supplementary measurement actions based on each component and the basic residual, connect the three-dimensional centerline of the underground pipeline according to the parameter uncertainty, and output the quality results. This is executed through the results evaluation module.

[0104] See attached document Figure 3 Step S100 includes steps S101, S102 and S103.

[0105] Step S101: Establish the correspondence between sensor spatial extrinsic parameters and unified sampling time. Rigid extrinsic parameter calibration was performed on the ground penetrating radar acquisition unit, the 3D point cloud acquisition unit, and the attitude acquisition unit to determine the fixed positional relationship between the phase center of the transmitting antenna, the phase center of the receiving antenna, the coordinate origin of the 3D point cloud acquisition unit, and the coordinate origin of the attitude acquisition unit.

[0106] After the extrinsic parameter calibration is completed, the 3D point cloud acquisition unit acquires point cloud identification markers on the rigid mounting frame. The distance between each point cloud identification marker is calculated based on its measured coordinates and compared with the pre-recorded calibration distance. If any distance difference exceeds the allowable deviation of the extrinsic parameter, the extrinsic parameter calibration is re-executed.

[0107] The time synchronization unit assigns a unified time reference to the ground-penetrating radar acquisition unit, the 3D point cloud acquisition unit, the attitude acquisition unit, the odometer acquisition unit, and the normal displacement unit. Each radar A-scan data record includes the radar sampling time, adjacent point cloud frame index, attitude data index, odometer data index, and normal displacement state index.

[0108] The initial time correspondence is used to establish a multi-source data index. The remaining time synchronization offset in forward and reverse measurements is not fixed in step S101, but is corrected in step S303 based on the directional odd component.

[0109] Step S102: Simultaneously acquire ground-penetrating radar data and 3D point cloud data along the area to be measured. The combined acquisition equipment moves along the road, hardened site, or soil surface to be measured. The ground-penetrating radar acquisition unit continuously acquires radar A-scan data according to the set sampling interval, while the three-dimensional point cloud acquisition unit simultaneously acquires three-dimensional point cloud data of the road and surrounding features.

[0110] The radar spatial sampling interval is determined based on the target's minimum pipe diameter, the radar center frequency, and the expected burial depth. Preferably, the radar spatial sampling interval is 0.005m to 0.100m. The average point spacing of the 3D point cloud near the radar trajectory is no greater than one-tenth of the radius of the fitted neighborhood on the ground surface.

[0111] The radar A-scan data undergoes DC component removal, early saturation suppression, background average response removal, out-of-band signal suppression, range gain compensation, and phase consistency processing. All processing parameters are kept consistent within the same local measurement area to avoid reflection arrival time shifts caused by parameter variations.

[0112] The coarse time-zero correction uses the direct wave initiation position or stable early response position as a reference to ensure that the early responses of different radar channels are within a unified time window. The coarse time-zero correction result is used for subsurface reflection candidate extraction and is not used as the final subsurface propagation time. The final subsurface propagation time is obtained through steps S302 and S303.

[0113] Step S103: Determine candidate areas for underground reflection based on the continuity of reflection time and spatial projection. Local energy extrema, phase extrema, continuous reflection ridges, and hyperbolic reflection patterns are extracted from radar B-scan data.

[0114] The reflection events in adjacent radar channels are matched. If the reflection time difference between two reflection events in adjacent radar channels is not greater than a preset trajectory time threshold, and the principal phase difference between the two reflection events is not greater than a preset trajectory phase threshold, the two reflection events are matched into the same reflection event sequence, and the matched reflection events form a continuous initial reflection trajectory along the survey line.

[0115] The preset trajectory time threshold is determined based on the radar sampling time resolution, the spacing between adjacent radar channels, and the expected maximum reflection time change rate. The preset trajectory phase threshold is determined based on the radar's effective bandwidth, signal-to-noise ratio, and phase extraction stability.

[0116] The trajectory of the joint acquisition device is determined based on attitude and odometer data. The initial reflection trajectory is then projected onto the ground surface represented by the 3D point cloud data according to the initial time index and the device trajectory.

[0117] For initial reflection trajectories in adjacent survey lines, when the distance between their surface projection positions is not greater than the preset survey line merging distance threshold and the reflection time difference is not greater than the preset survey line merging time threshold, the corresponding initial reflection trajectories are merged to form underground reflection candidate areas.

[0118] The width of the underground reflection candidate region is determined based on the discrete range of the reflection ridge line projected onto the surface and the initial positioning error. The length of the underground reflection candidate region is determined based on the distance between consecutive occurrences of reflection events. Isolated reflections with a continuous length below the minimum continuous length threshold are retained as block anomaly candidates and are not directly used as pipeline candidates.

[0119] An initial extension direction is fitted based on the center projection points of candidate underground reflection areas in adjacent survey lines. This initial extension direction is used in step S203 to deploy locally closed measurement units and is not used as the final pipeline routing output.

[0120] See attached document Figure 4 Step S200 includes steps S201, S202 and S203.

[0121] Step S201: Extract the surface point cloud of the candidate region and remove non-surface interference points. Based on the surface projection range of the candidate subsurface reflection region, the candidate region point cloud is extracted from the 3D point cloud data. The extracted range covers the candidate subsurface reflection region and the closed measurement baseline range on both sides.

[0122] Perform statistical outlier removal and radius outlier removal on the candidate region point cloud. Statistical outlier removal removes points whose deviation from the overall distribution exceeds a threshold based on the average distance from their neighbors. Radius outlier removal removes isolated points based on the number of neighboring points within a set spatial radius.

[0123] Surface point clouds are extracted based on the elevation continuity and local slope of the candidate region point clouds. Point clouds of vehicles, pedestrians, walls, guardrails, and green facilities are not included in the surface fitting. When the tops of manhole covers and valve boxes are continuous with the surrounding road surface, their top point clouds are included in the surface fitting, while their edge point clouds are recorded separately as boundary point clouds.

[0124] For point clouds obtained from forward and reverse measurements, spatial registration is performed using fixed curbstone, wall, manhole cover boundaries, and stable road surface features. If the average distance of the registered fixed feature point clouds exceeds the point cloud registration threshold, point cloud registration is re-performed. The preferred point cloud registration threshold is 0.005m to 0.050m.

[0125] Step S202: Calculate the local surface normal and establish local orthogonal measurement coordinates. In the For each surface sampling location, a neighborhood point cloud is selected. The mean coordinates and covariance matrix of the neighborhood point cloud are calculated. The unit eigenvector corresponding to the smallest eigenvalue of the covariance matrix is ​​then determined as the local surface normal. The formula for calculating the local surface normal is as follows: ; ; ; in, For the first The mean coordinates of the neighborhood point cloud corresponding to each surface sampling location represent the spatial geometric center of the local surface neighborhood, in meters. For the first The number of effective neighborhood points corresponding to each surface sampling location represents the number of point cloud points participating in local surface fitting, preferably 30 to 300. For the first The first surface sampling location corresponding to the first The coordinates of a neighboring point represent the three-dimensional position of the neighboring point in the engineering coordinate system, in meters. For the first The covariance matrix of the neighborhood point cloud corresponding to each surface sampling location represents the degree of dispersion of the neighborhood point cloud in three spatial directions, with units of m. 2 ; For the first Local surface normal unit vector corresponding to each surface sampling location; This is the minimum eigenvector extraction operator, which represents the extraction of the unit eigenvector corresponding to the minimum eigenvalue of the input matrix; This refers to the serial number of the surface sampling location; The neighboring point index, whose value ranges from 1 to 1. ; superscript This is the transpose notation for a matrix.

[0126] Adjust the local surface normal according to the positive elevation direction of the engineering coordinate system so that the local surface normal points upwards from the ground surface. When the angle between the local surface normal and the positive elevation direction is greater than 90°, reverse the local surface normal.

[0127] The initial extension direction of the candidate underground reflection region is fitted. In the local tangent plane perpendicular to the local surface normal, the projection of the initial extension direction onto the local tangent plane is normalized and determined as the longitudinal unit tangent. The unit vector perpendicular to the longitudinal unit tangent in the local tangent plane is determined as the transverse unit tangent.

[0128] The surface roughness is determined by the ratio of the minimum eigenvalue to the sum of all eigenvalues ​​of the neighborhood point cloud covariance matrix. If the surface roughness exceeds a roughness threshold, the neighborhood radius is expanded or reduced and the calculation is repeated. If the roughness threshold is exceeded in two consecutive calculations, the current surface sampling location is marked as an unsuitable normal displacement location.

[0129] Step S203: Deploy local closed measurement units according to the initial extension direction and the surface access boundary. The initial center projection position of the candidate underground reflection area is used as the center measuring point of the local closed measurement unit. Positive and negative transverse measuring points are respectively placed on both sides of the center measuring point along the transverse unit tangent. Positive and negative longitudinal measuring points are respectively placed on both sides of the center measuring point along the longitudinal unit tangent.

[0130] The design distances from the positive and negative lateral measuring points to the center measuring point are equal, as are the design distances from the positive and negative longitudinal measuring points to the center measuring point. The design lateral measurement baseline is determined based on the initial burial depth range, and the design longitudinal measurement baseline is determined based on the pipeline routing variation scale.

[0131] The closed-loop planning unit checks each measurement location based on the obstacle boundaries in the 3D point cloud data. When a measurement location falls into a vehicle, guardrail, green facility, curb face, or other impassable area, the corresponding symmetrical measurement points are moved synchronously to maintain the symmetry between the two measurement points about the central measurement point.

[0132] When the design symmetry cannot be maintained, record the adjusted actual measurement point coordinates and the actual measurement baseline. The closed-loop inversion module calculates the propagation geometric distance based on the actual measurement point coordinates, without performing nominal symmetry processing on the actual measurement point coordinates.

[0133] When the candidate area for underground reflection has two sets of intersecting extension directions, vertically separated reflection times, or obvious asymmetric reflection responses, four diagonal measuring points are added around the central measuring point. The diagonal measuring points are located at the positive transverse and positive longitudinal, positive transverse and negative longitudinal, negative transverse and positive longitudinal, and negative transverse and negative longitudinal positions, respectively.

[0134] See attached document Figure 5 Step S300 includes steps S301, S302 and S303.

[0135] Step S301: Conduct four-state observations of bidirectional and binormal distances at the same surface projection location. The joint acquisition module performs a forward first normal distance measurement at the measurement position determined in step S203. During the forward first normal distance measurement, the normal displacement unit keeps the phase center of the ground penetrating radar antenna at the first normal distance state.

[0136] The joint acquisition module performs a reverse first normal distance measurement along the opposite direction of movement. The reverse first normal distance measurement covers the same surface space as the forward first normal distance measurement. After completing the forward and reverse measurements of the first normal distance state, the normal displacement unit moves along the local surface normal to the second normal distance state.

[0137] The joint acquisition module sequentially performs forward and reverse second normal distance measurements under the second normal distance condition. All four observation conditions are completed within the same locally closed measurement unit.

[0138] Each state records the actual normal distance, attitude data, and odometer data, and determines the corresponding equipment trajectory based on the attitude data, odometer data, and 3D point cloud data. The ground-penetrating radar data is spatially resampled according to the 3D point cloud coordinates so that the four states correspond to the same surface projection location.

[0139] The distance difference between the first normal range state and the second normal range state is preferably 0.02m to 0.20m. If the radar reflection amplitude in the second normal range state is lower than the minimum signal-to-noise ratio requirement, the distance difference between the two normal range states is reduced and the data is re-acquired.

[0140] The time interval between completing the four observation states should preferably not exceed 30 minutes. If there is rainfall, road water accumulation, continuous ground vibration, or large-scale vehicle obstruction, the four-state observation should be suspended and resumed after the surface conditions have stabilized.

[0141] Step S302: Determine the air propagation path and correct the reflection arrival time based on the three-dimensional point cloud data. Obtain the extrinsic parameters between the pre-calibrated ground-penetrating radar antenna and the attitude sensor. Based on the 3D point cloud data, attitude data, extrinsic parameters, and local surface normal for each observation state, determine the phase center of the transmitting antenna and the phase center of the receiving antenna.

[0142] Along the local surface normal at the corresponding measurement location, the phase center of the transmitting antenna and the phase center of the receiving antenna are projected onto the ground surface to obtain the ground surface projection points of the transmitting end and the receiving end, respectively.

[0143] For a transceiver-based ground-penetrating radar, the sum of the distance from the phase center of the transmitting antenna to the surface projection point of the transmitting end and the distance from the phase center of the receiving antenna to the surface projection point of the receiving end is used as the air propagation path. For a transceiver-integrated ground-penetrating radar, twice the distance from the antenna phase center to the surface projection point is used as the air propagation path.

[0144] The formula for correcting the arrival time of underground reflections is as follows: ; in, For the first The candidate reflection event at the ... Measurement position and orientation status and normal distance state The corrected reflection arrival time after completing the air propagation time and the fixed delay correction inside the antenna, in ns; This represents the original reflection arrival time, in nanoseconds (ns). These are the coordinates of the phase center of the transmitting antenna, in meters (m). These are the coordinates of the transmitter's surface projection point, in meters. These are the coordinates of the phase center of the receiving antenna, in meters (m). These are the coordinates of the receiving end's surface projection point, in meters. Let be the speed of electromagnetic wave propagation in air, taken as 0.299792458 m / ns; The internal fixed delay of the pre-calibrated ground-penetrating radar antenna and acquisition circuit is expressed in nanoseconds (ns). This indicates the direction of motion, and its value is +1 or -1. This represents the normal distance state, with a value of 0 or 1. For measurement location sequence number; The candidate reflection event number; superscript Used to identify the transmitter; superscript Used to identify the receiving end; This represents the Euclidean distance in three-dimensional space.

[0145] After air propagation time correction, the ground-penetrating radar data under four observation conditions are resampled according to the three-dimensional point cloud coordinates and the surface projection position. The resampling interval is no greater than the joint initial radar spatial sampling interval. The reflection events under each condition are initially mapped according to the reflection principal phase, time window, and local waveform.

[0146] The internal fixed delay of the antenna is calibrated to the same value under the same operating conditions of the equipment. If the operating temperature of the equipment exceeds the calibrated temperature range or the antenna is replaced, the internal fixed delay of the antenna should be recalibrated.

[0147] Step S303: Perform orthogonal decomposition of the four-state reflection arrival time and execute trajectory synchronization correction and antenna coupling correction. For four corrected reflection arrival times at the same measurement location and for the same candidate reflection event, an orthogonal decomposition of the orientation state and normal distance state is performed. The four-state reciprocal decomposition formula is as follows: ; in, For the first The candidate reflection event at the ... The common reflection component at each measurement location, in nanoseconds; For directional odd components, the unit is ns; This represents the normal distance perturbation component, in nanoseconds (ns). The distance is the interaction component between the direction and the normal, and the unit is ns.

[0148] A synchronization correction model is constructed with the objective of minimizing the root mean square value of the directional odd components. The synchronization correction model uses the mileage ratio of reverse measurements, mileage translation, and time synchronization offset as the variables to be optimized. The synchronization correction model is solved, and the mileage ratio of reverse measurements, mileage translation, and time synchronization offset are updated based on the solution results. Spatial resampling and orthogonal decomposition of the four-state data are then performed again.

[0149] The termination condition for trajectory synchronization correction is that the root mean square value of the directional odd component of the valid candidate reflection event is not higher than the directional residual threshold, and the average spatial distance between the forward device trajectory and the reverse device trajectory is not higher than the point cloud registration threshold.

[0150] The direct wave or stable surface reference reflection is extracted from a preset time window starting at time zero. The system group delay is determined based on the time offset corresponding to the normal distance disturbance component of the reference reflection, and the system group delay is subtracted from the arrival time of each corrected reflection.

[0151] After antenna coupling correction is completed, the common reflection component, directional odd component, normal range perturbation component, and directional-normal range interaction component are recalculated and updated. The updated components overwrite the uncorrected components and serve as the sole input for subsequent steps.

[0152] The updated common reflection component is used as the stable reflection arrival time.

[0153] When the absolute value of the interaction component between the direction and normal distance exceeds the interaction threshold between the direction and normal distance, the corresponding reflection event is not directly included in the single pipeline common reflector observation group. The reflection event is included in the reflection set to be split, and is processed in step S403 in combination with the basic residual and diagonal measurement point data.

[0154] See attached document Figure 6 Step S400 includes steps S401, S402 and S403.

[0155] Step S401: Establish a common reflector observation group based on the arrival time, phase, and reciprocal components of stable reflections. The reflection correlation unit uses the common reflection component after completing trajectory synchronization correction and antenna coupling correction as the stable reflection arrival time.

[0156] Retrieve reflection events at different measurement locations. When the time difference between the stable reflection arrival times of two reflection events between adjacent measurement points is not greater than a preset continuity threshold, the two reflection events are determined to satisfy the time continuity condition.

[0157] The preset continuity threshold is determined based on the distance between measuring points, the expected pipeline burial depth, the range of local equivalent radar wave propagation speed, and the radar time resolution.

[0158] Calculate the principal phase difference between two reflection events that satisfy the time continuity condition. Extract the effective bandwidth from the spectrum of each reflection event and calculate the relative difference of the effective bandwidths between the two reflection events. The relative difference of the effective bandwidths is the ratio of the absolute value of the difference between the two effective bandwidths to the average value of the two effective bandwidths.

[0159] When the principal phase difference between two reflection events is not higher than the associated phase difference threshold, the relative difference in effective bandwidth is not higher than the bandwidth difference threshold, and the interaction component between direction and normal distance is not higher than the interaction threshold between direction and normal distance, an association relationship is established between the two reflection events.

[0160] The associated phase difference threshold is preferably 15° to 60°, the bandwidth difference threshold is preferably 0.10 to 0.50, and the interaction threshold between the direction and normal distance is preferably 0.03ns to 1.00ns.

[0161] A common reflector observation group is formed based on the correlation. Each common reflector observation group includes at least three valid reflection events corresponding to non-collinear measuring points. If there are fewer than three valid measuring points, the joint solution of pipeline position and local equivalent radar wave propagation velocity is not performed, and the process directly proceeds to step S502 to generate supplementary measurement actions.

[0162] The common reflector observation group records the stable reflection arrival time, the surface projection point at the transmitter, the surface projection point at the receiver, the main reflection phase, the effective reflection bandwidth, the directional odd component, the normal distance perturbation component, and the directional and normal distance interaction component for each reflection event.

[0163] In the initial association phase, the same reflection event can be included in multiple candidate common reflector observation groups. Step S402 calculates the candidate propagation velocity value and the fundamental residual for each candidate common reflector observation group, and step S403 determines the final assignment based on the closure result.

[0164] Step S402: Solve the pipeline parameters jointly based on the differential propagation time and symmetrical closure quantity. The first Each candidate pipe segment is represented as a spatial straight line, which is determined by a central control point and a directional unit vector. The central control point is the closest point from the axis of the candidate pipe segment to the center of a local closed measurement unit, and the central control point is located below the local ground surface.

[0165] For a separate transmit / receive ground-penetrating radar, the shortest distances from the surface projection points of the transmitter and receiver to the axis of the candidate pipe segment are calculated separately, and the sum of the two shortest distances is taken as the underground propagation geometric distance. For a transceiver integrated ground-penetrating radar, twice the shortest distance from the antenna surface projection point to the axis of the candidate pipe segment is taken as the underground propagation geometric distance.

[0166] The formula for calculating differential propagation time is as follows: ; ; ; in, For three-dimensional space points To the Candidate pipe section axis The shortest distance, in meters; It is a three-dimensional identity matrix; The direction unit vector of the candidate pipe segment; The coordinates of the center control point of the candidate pipe section are in meters. For the first The underground propagation geometric distance of each observation relative to the axis of the candidate pipe segment, in meters; For the first The first observation and the first The stable reflection arrival time difference between observations is expressed in nanoseconds. The local equivalent radar wave propagation speed is expressed in m / ns, preferably between 0.05 m / ns and 0.18 m / ns. This represents the residual error of differential propagation, expressed in nanoseconds (ns).

[0167] Observations within the same common reflector observation group share the same common system delay and the same pipeline radius term. By employing the stable reflection arrival time difference, the common system delay and pipeline radius term are eliminated from the differential propagation time relationship.

[0168] The transverse antisymmetric time difference and transverse second-order closure quantity are generated based on the stable reflection arrival times of the central measuring point, the positive transverse measuring point, and the negative transverse measuring point. The longitudinal antisymmetric time difference and longitudinal second-order closure quantity are generated based on the stable reflection arrival times of the central measuring point, the positive longitudinal measuring point, and the negative longitudinal measuring point.

[0169] The lateral antisymmetric time difference is the difference between the arrival time of stable reflection from the positive lateral measuring point and the arrival time of stable reflection from the negative lateral measuring point. The longitudinal antisymmetric time difference is the difference between the arrival time of stable reflection from the positive longitudinal measuring point and the arrival time of stable reflection from the negative longitudinal measuring point.

[0170] The second-order closure in the transverse direction is the sum of the stable reflection arrival times of the positive and negative transverse measuring points, minus twice the stable reflection arrival time of the center measuring point. The second-order closure in the longitudinal direction is generated in the same way from the positive, negative, and center measuring points.

[0171] Based on the differential propagation relationship between different measuring points, the antisymmetric time difference in the lateral and longitudinal directions, and the second-order closure quantities in the lateral and longitudinal directions, the local equivalent radar wave propagation velocity is calculated respectively, and the calculated values ​​are candidate values ​​for propagation velocity.

[0172] A joint objective function is established based on the differential propagation residual, the transverse antisymmetric closure residual, the longitudinal antisymmetric closure residual, the transverse second-order closure residual, the longitudinal second-order closure residual, and the interaction components of direction and normal distance. The joint solution and parameter covariance calculation formulas are as follows: ; ; in, For the first The vector of parameters to be determined for each candidate pipe segment; The optimal parameter vector; The differential observation weights range from 0 to 1. The residual for differential propagation is expressed in nanoseconds (ns). The residuals are lateral antisymmetric closed loops, expressed in nanoseconds (ns). The residuals are longitudinally antisymmetric closed loops, in nanoseconds (ns). The residuals are second-order closure errors in the transverse direction, expressed in nanoseconds (ns). The vertical second-order closed residual is expressed in nanoseconds (ns). The weighting coefficient for antisymmetric closed residuals is preferably between 0.1 and 20. The weighting coefficient for the second-order closed residual is preferably between 0.1 and 20. The weighting coefficient for the interaction components of direction and normal distance is preferably between 0.1 and 20. The parameter is the covariance matrix; This is the weighted residual variance, in nanoseconds (ns). 2 ; The residual Jacobian matrix; For the observation weight matrix; It is a regular matrix.

[0173] During the joint solution process, the length of the directional unit vector is constrained to 1, the central control point is constrained to be below the local ground surface, and the central control point is constrained to be the closest point from the candidate pipe segment axis to the center of the local closed measurement unit.

[0174] During the weighted robust optimization process, the observation weights are updated based on the residuals obtained in each iteration. The weights of observations that exceed the corresponding reference residual range are reduced until the parameters to be determined and the joint objective function satisfy the stopping condition.

[0175] The differential propagation residual, antisymmetric closure residual, and second-order closure residual are used to determine the basic residuals. The parameter uncertainty is determined by the parameter covariance matrix.

[0176] The initial value of the local equivalent radar wave propagation velocity is determined based on the radar reflection curvature or the propagation velocity in the reference area. After the joint solution is completed, the initial value is not used as the output value; instead, the output value is the optimal local equivalent radar wave propagation velocity corresponding to the joint objective function.

[0177] Between adjacent local pipe segments, continuous constraints are set for center position, continuous constraints for direction, and local smoothing constraints for propagation velocity. When the 3D point cloud data shows the boundaries of pavement material, backfill, or drainage structure, the local smoothing constraints for propagation velocity at the corresponding locations are reduced, so that the local equivalent radar wave propagation velocity is determined on both sides of the boundary.

[0178] Step S403: Confirm or split the common reflector observation group based on the candidate propagation velocity values, each component, and the fundamental residual. Based on the candidate pipe segment axis obtained in step S402, candidate propagation velocity values ​​are calculated using different effective observation pairs in the common reflector observation group. If the absolute value of the underground propagation geometric distance difference of the observation pair is lower than the minimum distance difference threshold, or the absolute value of the stable reflection arrival time difference is lower than the minimum time difference threshold, the observation pair is marked as a degraded observation pair and removed.

[0179] The mean and standard deviation of the propagation velocity candidate values ​​are calculated based on the valid observation pairs. When the differences in observation quality need to be considered, a weighted mean and weighted standard deviation are calculated using observation weights. The ratio of the standard deviation to the mean of the propagation velocity candidate values ​​is determined as the propagation velocity closure index.

[0180] Calculate the root mean square value of the interaction components of direction and normal distance within the common reflector observation group, and determine the root mean square value as the interaction closure index.

[0181] When the propagation velocity closure index is not higher than the propagation velocity closure threshold, the interaction closure index is not higher than the interaction threshold between direction and normal distance, and all basic residuals are not higher than their corresponding thresholds, the common reflector observation group is confirmed to be a single linear pipeline.

[0182] The propagation velocity closure threshold is determined based on the pipeline at a known location, the pipeline at the test site, or numerical simulation results, and is preferably between 0.03 and 0.20. A smaller value is used when the underground medium is homogeneous and the reflection signal-to-noise ratio is high, and a larger value is used when the underground medium is non-homogeneous and the pipeline is buried at a greater depth. The threshold for the interactive closure index is the interactive threshold between the direction and the normal distance.

[0183] When the propagation velocity closure index, the interaction closure index, or any fundamental residual exceeds the corresponding threshold, the target splitting unit divides the common reflector observation group into multiple subgroups based on the principal phase of the reflection, the arrival time of the stable reflection, the spatial projection position, and the interaction components of the direction and normal distance. Each subgroup then re-establishes the differential closure relationship and executes step S402.

[0184] When a locally closed measurement unit includes diagonal measurement points, the data from these diagonal points are used to verify whether there is mixed reflection between subgroups. If, after supplementary measurement and splitting, there are still no subgroups that satisfy the propagation velocity closure index, the interaction closure index, and the basic residual condition, the supplementary measurement action is generated again according to step S502.

[0185] If no subgroup meets the aforementioned conditions after reaching the preset maximum number of retests, the absolute value of the second-order closure quantity obtained by combining reflection events from different subgroups is determined as the mixed second-order closure residual. When there are diagonal measuring points, the second-order closure quantity obtained by combining reflection events corresponding to different diagonal measuring points and the central measuring point is preferred.

[0186] When the mixed second-order closure residual exceeds the diagonal closure threshold, the current common reflector observation group is marked as a nonlinear underground target group.

[0187] The diagonal closure threshold is determined by the statistical results of the mixed second-order closure residuals in the reference area without intersecting pipelines. It is preferably the sum of the mean of the mixed second-order closure residuals in the reference area and 2 to 4 times the standard deviation, and not less than 0.05 ns.

[0188] After target splitting, if two or more subgroups respectively satisfy the propagation velocity closure index, interaction closure index, and basic residual condition, the corresponding local pipe segments are output. When the center positions of different local pipe segments are close and their directions are basically parallel, they are marked as adjacent parallel pipelines. When the planar projections of different local pipe segments intersect and their elevations are different, they are marked as intersecting pipelines or pipelines stacked vertically.

[0189] See attached document Figure 7 Step S500 includes steps S501, S502 and S503.

[0190] Step S501: Calculate the residuals of each component, the fundamental residuals, the closure index, and the parameter uncertainty. The directional residual is calculated based on the directional odd components after secondary correction. The directional residual is the root mean square value of the directional odd components from valid observations within the same local candidate pipe segment. The directional residual is used to evaluate the consistency of mileage, time synchronization, and spatial trajectory between forward and reverse measurements.

[0191] The normal range residual is calculated based on the normal range perturbation component after secondary correction. The normal range residual uses the root mean square value of the normal range perturbation component from effective observations. The normal range residual is used to evaluate the remaining time difference between different normal range states after air propagation time correction and antenna coupling correction.

[0192] The cross-residuals of direction and normal distance are calculated based on the cross-components of direction and normal distance. These cross-residuals are derived from the root mean square values ​​of the cross-components of direction and normal distance from valid observations. The cross-residuals are used to evaluate whether changes in motion direction and normal distance can be explained separately.

[0193] The root mean square value is calculated based on the differential propagation residual from step S402, and the obtained root mean square value is determined as the differential propagation residual index.

[0194] Calculate the root mean square (RMS) value of the antisymmetric closure residual based on the horizontal and vertical antisymmetric closure residuals. Calculate the RMS value of the second-order closure residual based on the horizontal and vertical second-order closure residuals. The larger of the RMS value of the antisymmetric closure residual and the RMS value of the second-order closure residual is determined as the closure residual.

[0195] The differential propagation residual, antisymmetric closure residual, and second-order closure residual are determined as the basic residuals.

[0196] The standard deviations of the coordinate directions of the central control point, the pipeline burial depth, the pipeline orientation, and the local equivalent radar wave propagation velocity are extracted based on the parameter covariance matrix. The parameter uncertainty is characterized at least by the standard deviations of the planar position and burial depth.

[0197] When the standard deviation of pipeline burial depth is not higher than the standard deviation threshold of burial depth, and the standard deviation of planar position is not higher than the standard deviation threshold of planar position, the corresponding local pipeline section meets the surveying accuracy conditions.

[0198] The standard deviation thresholds for burial depth and planar position are set according to the accuracy requirements of the project. Preferably, the standard deviation threshold for burial depth is 0.03m to 0.20m, and the standard deviation threshold for planar position is 0.02m to 0.15m.

[0199] Step S502: Generate restricted retest actions based on different residual sources and re-acquire data. When the directional residual exceeds the directional residual threshold, the supplementary measurement planning unit maintains the current normal distance state and the current measurement baseline, and generates a reverse remeasurement action along the same trajectory. The joint acquisition module re-executes the reverse measurement along the original surface projection trajectory and shortens the time interval between the forward and reverse measurements.

[0200] When the normal distance residual exceeds the normal distance residual threshold, the supplementary measurement planning unit generates a supplementary measurement action to add a normal distance state. The added normal distance state can be a third normal distance state. The distance difference between the third normal distance state and the first normal distance state is different from the distance difference between the second normal distance state and the first normal distance state.

[0201] When adding a third normal distance state, auxiliary data acquisition is performed in the order of first normal distance state, third normal distance state, and first normal distance state. If the difference between the common reflection components of two first normal distance states exceeds the time drift threshold, device time drift correction is performed first, followed by normal distance response fitting.

[0202] The third normal distance state is used to determine the response relationship of the system group delay as a function of the normal distance. After the auxiliary acquisition is completed, two normal distance states that meet the signal-to-noise ratio and independence requirements are selected from the first, second, and third normal distance states, and bidirectional and bidirectional normal distance four-state observations are re-implemented in both the forward and reverse directions.

[0203] When the interaction residual between direction and normal distance exceeds the interaction threshold between direction and normal distance, the supplementary measurement planning unit generates a supplementary measurement action to add four diagonal measurement points and split the common reflector observation group. If necessary, the lateral or longitudinal measurement baseline is reduced. The supplementary measurement data enters step S303 to re-perform orthogonal decomposition, and then enters step S403 to perform the splitting of the common reflector observation group.

[0204] When the differential propagation residual index exceeds the differential propagation residual threshold, the supplementary measurement planning unit adjusts the center of the locally closed measurement unit according to the sign and magnitude of the differential propagation residual at the positive lateral, negative lateral, positive longitudinal, and negative longitudinal measurement points. The moving distance does not exceed half of the original lateral or longitudinal measurement baseline. Corresponding symmetrical measurement points are then re-established on both sides of the moved center measurement point.

[0205] When the closure residual exceeds the closure residual threshold, the supplementary measurement planning unit changes the lateral or longitudinal measurement baseline. When the antisymmetric closure residual component constituting the closure residual is large, the center of the local closure measurement unit is checked before changing the measurement baseline; when the second-order closure residual component is large, the measurement baseline is changed first. When the condition number of the parameter solution matrix exceeds the condition number threshold, the measurement baseline is also changed or a measurement line with a larger angle to the current pipeline direction is added.

[0206] The condition number threshold is preferably between 100 and 10000. When the condition number exceeds the threshold, it indicates that the current measurement geometry cannot adequately distinguish between pipeline burial depth and local equivalent radar wave propagation velocity.

[0207] After all the supplementary measurements were completed, the four-state observations of bidirectional and binormal distances were re-implemented, and the common reflector observation group and differential closure relationship were re-established.

[0208] The preset maximum number of retests for a single local closure measurement unit is preferably 2 to 6. If the closure condition and accuracy condition are still not met after reaching the preset maximum number of retests, an unclosed anomaly mark will be output for the current area.

[0209] Step S503: Connect the axes of adjacent candidate pipe segments and output the three-dimensional centerline and quality results of the underground pipeline. Adjacent candidate pipe segment axes obtained along the same underground reflection candidate region are connected. Each candidate pipe segment axis has two endpoints within the corresponding range of a locally closed measurement unit. The distance between two endpoints that are close to each other in adjacent candidate pipe segment axes is the adjacent endpoint distance.

[0210] The connection conditions include at least the following: the distance between adjacent endpoints is not greater than a preset connection distance threshold, the directional angle is not greater than a preset connection angle threshold, each closure index does not exceed its corresponding threshold, each residual does not exceed its corresponding threshold, and the standard deviation of the plane position and the standard deviation of the burial depth do not exceed their corresponding thresholds.

[0211] Furthermore, it can be required that the difference in burial depth between adjacent candidate pipe sections does not exceed the maximum difference in burial depth, and the difference in local equivalent radar wave propagation velocity does not exceed the maximum velocity difference.

[0212] The connection distance threshold is determined based on the spacing between local closed measurement units, preferably 1 to 2 times the spacing between local closed measurement units. The connection angle threshold is preferably 5° to 35°. The maximum burial depth difference is determined based on the underground pipeline slope requirements and the spacing between adjacent local closed measurement units.

[0213] When all closure indices, residuals, and standard deviations do not exceed their corresponding thresholds, and the distance between adjacent endpoints and the angle between their directions are not greater than the connection distance threshold and the connection angle threshold, respectively, the axes of adjacent candidate pipe segments are connected to form the three-dimensional centerline of the underground pipeline.

[0214] When manholes, valve boxes, or well chamber boundaries are identified in 3D point cloud data, these boundaries are used as candidate regions for topology connections. When a local pipe segment extends into a candidate region for topology connections and its radar reflection is continuous in the vicinity of the candidate region, a topology connection is established between the local pipe segment and the corresponding well chamber.

[0215] When adjacent local pipe segments intersect on the plane projection but have different elevations, and each satisfies the propagation velocity closure condition, the three-dimensional centerlines of the two underground pipelines remain independent, and no node merging is performed. When the centers of adjacent local pipe segments are close in position and have the same direction, but the reflected principal phase and the local equivalent radar wave propagation velocity are different, they are maintained as different pipeline candidate results.

[0216] The output unit outputs the three-dimensional centerline, coordinates of the center control point, local direction, pipeline burial depth, pipeline center elevation, local equivalent radar wave propagation velocity, standard deviation of plane position, standard deviation of burial depth, propagation velocity closure index, and original radar observation index for each underground pipeline.

[0217] For local pipe sections that meet all closure and accuracy requirements, valid mapping results are output. For local pipe sections that only partially meet the requirements, mapping results with quality limitations are output. For areas that still do not meet the single-pipe closure conditions after reaching the preset maximum number of re-measurements, the following conditions are output: multiple pipes not separated, strong clutter, unstable antenna coupling, unidentifiable propagation speed, or requiring verification by other detection methods. An unclosed anomaly marker is also output.

[0218] Specific application examples: See attached document Figure 8 , Figure 8 The horizontal axis represents the iteration number, ranging from 0 to 6; the left vertical axis represents the common reflection time, in nanoseconds (ns); and the right vertical axis represents the residual, also in nanoseconds (ns). A dashed grid is used in the figure to aid in reading the curve values.

[0219] Figure 8 In the diagram, solid dotted curves represent common reflection components, dashed square dotted curves represent directional residuals, dotted-dash triangular dotted curves represent the interaction residuals between directional and normal distances, and dotted-line diamond dotted curves represent differential propagation residual indices. Figure 8The system also sets directional residual threshold lines and interactive residual threshold lines. The directional residual threshold line is represented by a long gray dashed line with a threshold of 0.20 ns; the interactive residual threshold line is represented by a gray dotted line with a threshold of 0.08 ns.

[0220] This embodiment selects a 180m long road with an average width of 12m in a county's urban area as the measurement area. The road surface is asphalt pavement, with manholes, valve boxes, and road repair strips in some areas. The measurement area includes a known control section and an unknown detection section. The known control section is used to determine the antenna's internal fixed delay, residual threshold, and mapping result deviation, while the unknown detection section is used to perform underground pipeline mapping and identification.

[0221] This embodiment uses a ground-penetrating radar with a center frequency of 400MHz. The radar spatial sampling interval is set to 0.02m, and the average point spacing of the three-dimensional point cloud acquisition unit near the radar trajectory is 0.015m. The first normal distance state is set to 0.06m, the second normal distance state is set to 0.12m, and the distance difference between the two normal distance states is 0.06m.

[0222] Joint preliminary surveys were conducted along the longitudinal and transverse directions of the road. The joint preliminary surveys identified three continuous candidate areas for underground reflection. The first candidate area for underground reflection extends longitudinally along the road; the second candidate area for underground reflection intersects with the planar projection of the first candidate area for underground reflection near the manhole; the third candidate area for underground reflection is locally concentrated and does not form a continuous linear reflection trajectory.

[0223] A locally closed measurement unit was set up in the first candidate area for underground reflection. The lateral measurement baseline was set to 0.60m, and the longitudinal measurement baseline was set to 0.80m. The neighborhood radius of the point cloud of the center measuring point, positive lateral measuring point, negative lateral measuring point, positive longitudinal measuring point, and negative longitudinal measuring point was all set to 0.35m.

[0224] The number of effective neighboring points obtained around the central measuring point is: Coordinates of 120 neighboring points Substituting into step S202, the local surface normal calculation process is used to obtain the neighborhood point cloud covariance matrix. The three eigenvalues ​​are 0.000012m. 2 0.0064m 2 and 0.0091m 2 .

[0225] covariance matrix The local surface normal unit vector corresponding to the minimum eigenvalue for: ; Local surface normal unit vector The elevation direction of the ground-penetrating radar antenna does not completely coincide with the engineering coordinate system elevation direction. When directly raising and lowering the antenna using the engineering coordinate system elevation direction, the surface projection point corresponding to the second normal distance state is offset by 0.004m relative to the first normal distance state. A local surface normal unit vector is used. After controlling the normal displacement unit, the surface projection point corresponding to the two normal distance states shifts by 0.001m. Therefore, the horizontal projection deviation in the dual normal distance observation is reduced by 0.003m.

[0226] In the positive first normal distance observation at the central measuring point, the first The original reflection arrival time of each candidate reflection event The phase center of the transmitting antenna is 21.68 ns. Projection point to the ground surface of the launch end The distance is 0.12m, and the phase center of the receiving antenna is... To the ground projection point of the receiving end The distance is 0.12m.

[0227] Speed ​​of electromagnetic waves in air Take 0.299792458 m / ns, and use the fixed internal delay of the antenna. Take 0.06 ns. Airborne propagation time. The calculation results are as follows: ; Arrival time of the original reflection Airborne transmission time and fixed delay inside the antenna Substituting the correction relationship from step S302, we obtain the corrected reflection arrival time under the positive first normal distance condition: ; Based on the sampling time resolution of ground penetrating radar, the corrected reflection arrival time will be... The time was recorded as 20.82 ns.

[0228] The corrected reflection arrival times corresponding to the reverse first normal distance state, the forward second normal distance state, and the reverse second normal distance state at the same central measuring point are as follows: ; ; ; Substitute the arrival times of the four corrected reflections into the four-state reciprocal decomposition process in step S303. Common reflection component. The calculation results are as follows: ; directional odd components The calculation results are as follows: ; Normal distance perturbation component The calculation results are as follows: ; Interaction components of direction and normal distance The calculation results are as follows: ; In this embodiment, based on the statistical results of the known control segment reference area, the directional residual threshold is set to 0.20 ns, the normal distance residual threshold is set to 0.15 ns, and the directional and normal distance interaction threshold is set to 0.08 ns.

[0229] The directional odd component of the center measuring point Below the directional residual threshold of 0.20 ns; normal range perturbation component The residual value is below the normal distance residual threshold of 0.15 ns; the interaction component between direction and normal distance. The distance between the direction and normal is below the interaction threshold of 0.08 ns. Therefore, candidate reflection events are retained as valid reflection events, and the common reflection component is... This serves as the stable reflection arrival time at the measuring point.

[0230] Before trajectory synchronization correction, the directional residual of all valid measuring points in the first underground reflection candidate area was 0.42 ns, which was higher than the directional residual threshold of 0.20 ns. The reciprocity correction module solved the synchronization correction model based on the directional odd components, adjusted the mileage ratio, mileage translation, and time synchronization offset of the reverse measurement, and resampled the forward and reverse measurement data spatially.

[0231] Figure 8 In the diagram, an iteration count of 0 indicates the state of four-state decomposition using the initial mileage relationship and the initial time synchronization relationship; iteration counts of 1 to 6 indicate the states after completing the first to sixth trajectory synchronization corrections and adjusting the closed inversion parameters, respectively.

[0232] Figure 8 The values ​​of the common reflection components shown are as follows for iterations 0 to 6: ; The common reflection component gradually increased from the initial 20.61 ns, reaching 20.75 ns in the fourth iteration, and remained at 20.75 ns in the fifth and sixth iterations. The fact that the common reflection component maintained the same result for three consecutive iterations indicates that the effects of air propagation time correction, trajectory synchronization correction, and closure inversion parameter adjustment on the stable reflection arrival time have converged.

[0233] Figure 8 The values ​​of the directional residuals shown are as follows for iterations 0 to 6: ; The directional residual was 0.23 ns in the second iteration, still above the 0.20 ns threshold; it decreased to 0.17 ns in the third iteration, falling below the threshold for the first time. The reciprocity correction module continued to perform small-scale mileage ratio and time synchronization offset adjustments, further reducing the directional residual to 0.145 ns, 0.136 ns, and 0.130 ns in the fourth to sixth iterations, respectively.

[0234] Therefore, the third iteration was used to confirm that the directional residual met the threshold condition, and the sixth iteration corresponds to the final stable result. The final directional residual was 0.13 ns, a decrease of 0.29 ns compared to the initial directional residual of 0.42 ns, representing a reduction of 69.05%.

[0235] Figure 8 The values ​​of the interaction residuals between the direction and normal distance shown in the figure for iterations 0 to 6 are as follows: ; The interaction residual between the direction and normal distance was 0.19 ns in the initial state, which was higher than the interaction threshold of 0.08 ns. In the first and second iterations, it decreased to 0.15 ns and 0.11 ns, respectively, but was still higher than the interaction threshold. In the third iteration, it decreased to 0.075 ns, which was lower than the interaction threshold of 0.08 ns for the first time.

[0236] In the 4th to 6th iterations, the interaction residuals between the direction and normal distances were further reduced to 0.058 ns, 0.046 ns, and 0.040 ns, respectively. The final interaction residuals between the direction and normal distances decreased by 0.15 ns relative to the initial value, representing a reduction of 78.95%, indicating that the non-independent response between the forward and reverse measurement differences and the change in normal distance was suppressed.

[0237] Figure 8 The values ​​of the differential propagation residual index shown are as follows for iterations 0 to 6: ; The differential propagation residual index continuously decreases with the joint updates of the candidate pipe segment center position, pipeline direction, and local equivalent radar wave propagation velocity. In the 4th iteration, the differential propagation residual index decreases to 0.22 ns; in the 5th iteration, it decreases to 0.19 ns; and in the 6th iteration, it decreases to 0.18 ns. The change from the 5th to the 6th iteration is 0.01 ns, which is lower than the iteration change threshold set in this embodiment, therefore the joint iteration is stopped.

[0238] Figure 8The directional residual threshold line remained at 0.20 ns, and the interaction residual threshold line remained at 0.08 ns. Both the directional residual curve and the interaction residual curve between the direction and normal distance fell below their respective threshold lines after the third iteration. The common reflection component curve remained horizontal after the fourth iteration, and the differential propagation residual index curve tended to stabilize after the fifth iteration. These curves collectively constitute the termination criteria for trajectory synchronization correction and closure inversion.

[0239] The common reflector observation group for the first underground reflection candidate area includes the effective reflection events corresponding to the central measuring point, positive transverse measuring point, negative transverse measuring point, positive longitudinal measuring point, and negative longitudinal measuring point.

[0240] Among them, the One valid observation relative to the candidate pipe segment axis underground propagation geometric distance It is 2.10m, the first One valid observation relative to the candidate pipe segment axis underground propagation geometric distance It is 2.34m.

[0241] No. The first valid observation and the first Stable reflection arrival time difference between valid observations It is -2.40 ns. Based on the differential propagation time relationship in step S402, the locally equivalent radar wave propagation velocity is... for: ; The remaining valid observation pairs in the common reflector observation group were calculated in the same way, and the obtained candidate propagation velocity values ​​ranged from 0.094 m / ns to 0.106 m / ns.

[0242] The closed inversion module propagates the residuals from the differences of each valid observation pair. Lateral antisymmetric closed residual Longitudinal antisymmetric closed residual Lateral second-order closed residual Longitudinal second-order closed residual Interaction components of direction and normal distance Input the joint objective function.

[0243] This embodiment uses the antisymmetric closed residual weighting coefficient Setting it to 4 will adjust the weighting coefficient of the second-order closed residual. Setting it to 6 will adjust the weighting coefficients of the direction and normal distance interaction components. Set it to 3.

[0244] The initial pipeline location, initial pipeline direction, and local equivalent radar wave propagation velocity of 0.10 m / ns obtained from the joint preliminary survey are used as the parameter vector to be determined. The initial values ​​are determined. After 6 iterations, the relative change in the joint objective function falls below the stopping threshold, yielding the optimal parameter vector. .

[0245] Optimal parameter vector The corresponding local pipe segment center control point has a lateral offset of 0.08m and a longitudinal offset of -0.03m relative to the center of the local closed measurement unit. The normal burial depth is 1.05m, and the local equivalent radar wave propagation velocity is... It is 0.099 m / ns.

[0246] Calculate the parameter covariance matrix based on the parameter covariance relationship in step S402. From the parametric covariance matrix The standard deviation of the obtained planar position is 0.04m, the standard deviation of the burial depth is 0.05m, and the standard deviation of the local equivalent radar wave propagation velocity is 0.006m / ns.

[0247] In this embodiment, the standard deviation threshold for planar position is set to 0.08m, the standard deviation threshold for burial depth is set to 0.10m, and the standard deviation threshold for local equivalent radar wave propagation velocity is set to 0.012m / ns. Since all three standard deviations for the first underground reflection candidate area are below the corresponding thresholds, the local pipe segment mapping results for the first underground reflection candidate area are output.

[0248] The initial common reflector observation set for the second candidate underground reflection region includes multiple reflection events near the inspection well. The initially calculated interaction closure index is 0.19 ns, which is higher than the 0.08 ns interaction threshold between direction and normal distance.

[0249] Based on the residual type corresponding to the interactive closure index, the result evaluation module does not perform reverse retesting of the same trajectory, nor does it directly output a single pipeline. Instead, it generates a supplementary measurement action that adds 4 diagonal measuring points and splits the common reflector observation group, while adjusting the lateral measurement baseline from 0.60m to 0.40m.

[0250] After completing the diagonal measurement point supplementary measurement, the reflection correlation unit divides the initial common reflector observation group into two subgroups based on the main phase of the reflection, the arrival time of the stable reflection, the spatial projection position, and the interaction components of the direction and normal distance.

[0251] The mean local equivalent radar wave propagation velocity of the first subgroup is 0.102 m / ns, and the propagation velocity closure index is 0.07. The mean local equivalent radar wave propagation velocity of the second subgroup is 0.087 m / ns, and the propagation velocity closure index is 0.09.

[0252] In this embodiment, the propagation velocity closure threshold is set to 0.15. The propagation velocity closure indices of both the first and second subgroups are below 0.15, and the interaction closure indices after the two subgroups have split are reduced to 0.05 ns and 0.04 ns, respectively, both below the interaction threshold of 0.08 ns for directional and normal distance. The basic residuals of both subgroups are also below the corresponding thresholds.

[0253] The pipeline planar projections of the first and second subgroups intersect near the inspection well, with a center elevation difference of 0.43m between the two subgroups. Since the two subgroups satisfy the propagation velocity closure index, the interaction closure index, and the basic residual condition respectively, and their center elevations are different, the closure inversion module recognizes the two subgroups as two independent underground pipelines and does not merge the centerline nodes at the intersection.

[0254] The third candidate area for underground reflection showed strong, localized blocky reflections in the joint preliminary measurement. Since the initial common reflector observation group could not meet the propagation velocity closure condition, the results evaluation module generated a supplementary measurement action at the diagonal measurement points.

[0255] After completing the diagonal measurement point supplementary measurements, the mixed second-order closure residual was calculated based on the second-order closure values ​​obtained from different subgroup reflection event combinations. The mixed second-order closure residual for the third underground reflection candidate region was 0.61 ns, which is higher than the diagonal closure threshold of 0.20 ns. The propagation velocity candidate values ​​calculated from different effective observation pairs ranged from 0.052 m / ns to 0.168 m / ns, with a propagation velocity closure index of 0.38, which is higher than the propagation velocity closure threshold of 0.15.

[0256] After the second supplementary measurement of the third underground reflection candidate area, there was still no subgroup that simultaneously met the propagation velocity closure index, the interaction closure index, and the basic residual condition, and the preset maximum number of supplementary measurements was reached. Its mixed second-order closure residual was still higher than the diagonal closure threshold, and the propagation velocity closure index was still higher than the propagation velocity closure threshold. The output unit marked the third underground reflection candidate area as a non-linear underground target group and did not connect it to the three-dimensional centerline of the underground pipelines in the first or second underground reflection candidate areas.

[0257] Appendix Figure 8 In the process, the common reflection component curve remained at 20.75 ns after the fourth iteration, the direction residual curve was below 0.20 ns after the third iteration, the cross-residual curve of direction and normal distance was below 0.08 ns after the third iteration, and the differential propagation residual index curve decreased from 0.74 ns to 0.18 ns. Therefore, Figure 8 It fully illustrates the process of stable reflection arrival time, reduced trajectory synchronization error, reduced interaction error between direction and normal distance, and gradual convergence of the differential propagation model.

[0258] Experimental verification and effect comparison: See attached document Figure 9 , Figure 9 The horizontal axis represents the verification point number, ranging from 1 to 30; the vertical axis represents the absolute error of burial depth, in meters. Figure 9 A dashed background grid was used to assist in reading the error of each verification point. The area where verification points 10 and 11 are located is represented by the first gray area and marked as a road repair strip; the area where verification points 24 and 25 are located is represented by the second gray area and marked as a pipeline intersection area.

[0259] Figure 9 In the diagram, the solid dot curve represents the fixed propagation speed method, the dashed square dot curve represents the point cloud coordinate fusion method, the dotted-dash triangular dot curve represents the single normal distance forward and reverse correction method, and the dotted-line diamond dot curve represents the method of this invention.

[0260] The experimental verification area was set up within the scope of a road renovation project in a certain district. The verification area included 30 excavation verification points, of which 18 were for metal pipelines and 12 were for non-metal pipelines. The pipeline burial depth ranged from 0.55m to 2.10m, and the pipe diameter ranged from 0.08m to 0.80m.

[0261] The 30 verification points include 12 single-pipe locations, 8 parallel-pipe locations, 6 intersecting-pipe locations, and 4 locations with strong local clutter. After excavation, engineering surveying equipment was used to obtain the horizontal position and elevation of the pipeline center, which were then used as verification benchmarks.

[0262] To reduce the impact of vehicle vibration, temporary obstruction, and equipment start-up drift on the results during a single data collection, each verification point was collected three times using the same equipment parameters. Figure 9 For each curve point, the median value of the absolute error of the burial depth is taken from the three measurement results of the corresponding verification point. The root mean square error of the plane position, the root mean square error of the burial depth, and the maximum error are calculated based on the original results of all repeated measurements.

[0263] The experiment included four sets of processing results: a fixed propagation velocity method, a point cloud coordinate fusion method, a single normal distance forward and reverse correction method, and the method of this invention. All four methods used the same raw ground-penetrating radar data, the same 3D point cloud data, and the same candidate subsurface reflection regions.

[0264] The fixed propagation velocity method uses a uniform propagation velocity of 0.10 m / ns and calculates the pipeline burial depth based on the arrival time of reflection from a single survey line. The point cloud coordinate fusion method uses three-dimensional point clouds to correct radar trajectories and surface elevations, but does not perform dual-normal range acquisition or differential closure inversion.

[0265] The single normal range forward and reverse correction method performs forward and reverse measurements and corrects the mileage and time synchronization offset based on the difference between the forward and reverse measurements, but does not change the antenna normal range and does not calculate the normal range disturbance component and the interaction component between the direction and the normal range.

[0266] The method of the present invention is executed according to steps S100 to S500. The first normal distance state is 0.05m to 0.08m, the second normal distance state is 0.10m to 0.16m, the lateral measurement baseline is 0.40m to 0.90m, and the longitudinal measurement baseline is 0.50m to 1.20m.

[0267] Figure 9 The median absolute error of burial depth at verification points 1 to 30 using the fixed propagation velocity method are as follows: ; The median absolute error of burial depth for the fixed propagation velocity method ranges from 0.14m to 0.33m. The errors at verification points 4, 10, 17, 20, 25, and 28 are 0.31m, 0.33m, 0.29m, 0.30m, 0.31m, and 0.28m, respectively, forming several distinct peaks.

[0268] Figure 9 The median absolute error of burial depth at verification points 1 to 30 using the midpoint cloud coordinate fusion method are as follows: ; The median absolute error of the burial depth for the point cloud coordinate fusion method ranges from 0.07m to 0.20m. Compared with the fixed propagation speed method, the overall error of each verification point is reduced, but the errors at verification points 10 and 25 still reach 0.20m and 0.19m, respectively.

[0269] Figure 9 The median values ​​of the absolute burial depth error for the single normal distance forward and reverse correction method at verification points 1 to 30 are as follows: ; The median absolute error of the burial depth for the single-normal distance forward and reverse correction method is distributed between 0.05m and 0.15m. The error at verification points 10 and 25 is 0.15m, indicating that forward and reverse correction can reduce mileage and synchronization errors, but it cannot independently correct for antenna coupling changes and multi-target mixed responses in road repair strips and pipeline intersection areas.

[0270] Figure 9 The median absolute error of burial depth at verification points 1 to 30 of the present invention is as follows: ; The median absolute error of the burial depth illustrated in the figure for the method of the present invention ranges from 0.03m to 0.09m. The curves for the 30 verification points are all below the corresponding curves for the other three methods. The errors at verification points 10 and 25 are both 0.09m, which are the larger values ​​in the figure for the method of the present invention, but are still 0.24m and 0.22m lower than those for the fixed propagation speed method, respectively.

[0271] The road repair strip corresponds to verification points numbered 10 and 11. Within this area, the median absolute error of burial depth for the fixed propagation velocity method is 0.33m and 0.17m, respectively; for the point cloud coordinate fusion method, it is 0.20m and 0.09m, respectively; for the single normal distance forward and reverse correction method, it is 0.15m and 0.07m, respectively; and for the method of this invention, it is 0.09m and 0.04m, respectively.

[0272] Verification point number 10 represents the location of road repair material changes, where the error using the fixed propagation velocity method reached 0.33m. After point cloud coordinate fusion, this error was reduced to 0.20m; after single-normal distance forward and reverse correction, it was reduced to 0.15m; and after using the method of this invention, it was reduced to 0.09m. This result demonstrates that the joint solution of local surface normal correction and local equivalent radar wave propagation velocity can suppress continuous burial depth deviations caused by changes in road repair materials.

[0273] The pipeline crossing area corresponds to verification points 24 and 25. Within this area, the median absolute error of burial depth using the fixed propagation velocity method is 0.22m and 0.31m, respectively; the median error using the point cloud coordinate fusion method is 0.12m and 0.19m, respectively; the median error using the single normal distance forward and reverse correction method is 0.09m and 0.15m, respectively; and the median error using the method of this invention is 0.05m and 0.09m, respectively.

[0274] Verification point number 25 is the location where the planar projections of two underground pipelines intersect. The error of the fixed propagation velocity method is 0.31m, the error of the point cloud coordinate fusion method is 0.19m, the error of the single normal distance forward and reverse correction method is 0.15m, and the error of the method of this invention is 0.09m. This result shows that the diagonal measurement point closure relationship and the splitting of the common reflector observation group can reduce the positioning error caused by the mixing of reflections from intersecting pipelines.

[0275] The root mean square error of the planar position and the root mean square error of the burial depth were 0.152 m and 0.226 m respectively, based on all repeated validation results using the fixed propagation velocity method. The burial depth error curves for road repair strips and areas with varying moisture content showed continuous shifts.

[0276] The root mean square error of the planar position and the root mean square error of the burial depth using the point cloud coordinate fusion method are 0.108m and 0.154m, respectively. This method reduces the elevation error caused by surface undulations, but burial depth deviations still exist in areas where the propagation velocity varies.

[0277] The root mean square error of the planar position and the root mean square error of the burial depth using the single normal range forward and reverse correction method are 0.086 m and 0.118 m, respectively. This method reduces the planar position error caused by the slippage of the ranging wheel and the time synchronization offset, but reflection time offset still exists in the area where the antenna's altitude changes.

[0278] The root mean square error of the planar position of the method of this invention is 0.061m, and the root mean square error of the burial depth is 0.073m. The maximum planar position error among all repeated measurements is 0.128m, and the maximum burial depth error is 0.146m. Figure 9 The curve points shown are the median values ​​of three measurements at each verification location. Therefore, the maximum value of a single point and the maximum single error in all repeated measurements are statistically analyzed separately.

[0279] Compared with the fixed propagation speed method, the root mean square error of the planar position of the method of the present invention is reduced by 0.091m, a reduction of 59.87%; and the root mean square error of the burial depth is reduced by 0.153m, a reduction of 67.70%.

[0280] Compared with the point cloud coordinate fusion method, the root mean square error of the planar position of the method of the present invention is reduced by 0.047m, which is 43.52%; the root mean square error of the burial depth is reduced by 0.081m, which is 52.60%.

[0281] Compared with the single normal distance forward and reverse correction method, the root mean square error of the planar position of the method of the present invention is reduced by 0.025m, a reduction of 29.07%; the root mean square error of the burial depth is reduced by 0.045m, a reduction of 38.14%.

[0282] In the eight parallel pipeline verification locations, the fixed propagation velocity method merges the two pipelines at three locations into a single reflective target, the point cloud coordinate fusion method merges the two pipelines at two locations into a single reflective target, and the single normal distance forward and reverse correction method merges the two pipelines at one location into a single reflective target. The method of this invention uses propagation velocity closure and interactive closure indices for splitting, outputting two independent three-dimensional centerlines for each of the eight locations.

[0283] In the six cross-pipeline verification locations, the method of this invention forms two common reflector observation groups for five locations, satisfying the propagation velocity closure index, the interaction closure index, and the basic residual condition. For the remaining location, because the elevation difference between the pipeline centers is lower than the radar's time resolution capability, after reaching the preset maximum number of supplementary measurements, the output shows a state of multiple pipelines not being separated, and no single, definitive pipeline result is output.

[0284] In the four verification locations with strong local clutter, the fixed propagation velocity method interprets two of the blocky reflections as short pipelines, while the point cloud coordinate fusion method interprets one of the blocky reflections as a short pipeline. In the method of this invention, the interactive closure index or mixed second-order closure residuals at all four locations exceed the corresponding thresholds, and all outputs are the state of the nonlinear underground target group.

[0285] Figure 9 The four curves all show local increases near verification points 4, 10, 14, 17, 20, 25, and 28, indicating that these locations correspond to relatively complex surface conditions, local medium variations, or multi-target reflection conditions. The curve trends of the four methods are basically consistent, but the dotted-line diamond curve of the present invention is always below the other curves, indicating that the present invention can maintain a low absolute error in burial depth at the same complex locations.

[0286] The curve of the fixed propagation speed method exhibits the largest fluctuation range, while the curve fluctuations of the point cloud coordinate fusion method and the single normal distance forward and reverse correction method decrease sequentially. The curve of the method of this invention has the smallest value range, and the road repair strip and pipeline intersection area do not form the same high amplitude peak value as the fixed propagation speed method.

[0287] The workload of supplementary measurements at 30 verification points was statistically analyzed. The total length of the supplementary measurement line using the regular grid supplementary measurement method was 486m. The method of this invention performs directional supplementary measurement based on the residual category, reducing the total length of the supplementary measurement line to 297m, a reduction of 189m, or 38.89%.

[0288] In the supplementary measurement actions generated by the method of the present invention, reverse re-measurement actions on the same trajectory account for 31% of all supplementary measurement actions, actions to increase the normal distance state account for 18%, actions to adjust the center of the local closed measurement unit account for 24%, actions to change the measurement baseline account for 17%, and actions to add diagonal measurement points and split the common reflector observation group account for 10%.

[0289] The directional residuals corresponding to the reverse retest of the same trajectory decreased by an average of 0.21 ns after the retest. The normal distance residuals corresponding to the action of increasing the normal distance state decreased by an average of 0.17 ns after the retest. The directional and normal distance interaction residuals corresponding to the action of adding diagonal measuring points and splitting the common reflector observation group decreased by an average of 0.23 ns after the target split was completed.

[0290] Experimental results show that after the local surface normal determined by the 3D point cloud data is used to constrain the antenna normal displacement and air propagation path, the reflection time correction results of the slope and road repair areas remain continuous. The bidirectional and binormal distance four-state observation outputs the trajectory synchronization response, antenna coupling response, and directional and normal distance interactive response separately, allowing different error sources to correspond to different remeasurement actions.

[0291] Appendix Figure 8 The results show that the common reflection component stabilizes at 20.75 ns with iteration, the directional residual decreases from 0.42 ns to 0.13 ns, the cross-residual between directional and normal distances decreases from 0.19 ns to 0.04 ns, and the differential propagation residual index decreases from 0.74 ns to 0.18 ns. Both the directional residual and the cross-residual between directional and normal distances fall below their respective thresholds after the third iteration.

[0292] Appendix Figure 9 The results show that the error variation trends of the fixed propagation velocity method, the point cloud coordinate fusion method, the single normal distance forward and reverse correction method, and the method of the present invention are correlated at the same verification point. However, the median absolute error of the burial depth of the present invention at 30 verification points is lower than that of the other three methods. The error reduction results in road repair strips and pipeline intersection areas indicate that the method of the present invention can simultaneously suppress the positioning deviation caused by local propagation velocity changes, antenna coupling changes, and multi-target reflection mixing.

[0293] Differential propagation time and symmetric closure eliminate common system delay and pipeline radius terms within the same locally closed measurement unit. The joint solution determines the locally equivalent radar wave propagation velocity based on local propagation conditions. The propagation velocity closure index, interactive closure index, and fundamental residual are used together to constrain the common reflector observation group. Observations that do not meet the closure conditions are not output as determining the three-dimensional centerline of the pipeline.

Claims

1. A method for mapping and identifying underground pipelines using ground-penetrating radar combined with 3D point cloud, characterized in that, Includes the following steps: Acquire ground-penetrating radar data, 3D point cloud data, attitude data, and odometer data, and determine candidate areas for underground reflection; Calculate the local surface normal based on the 3D point cloud data within the candidate area and deploy local closed measurement units; Four-state observations of bidirectional and binormal distances were implemented in the measurement unit. Correct the reflection arrival time and orthogonally decompose to obtain the common reflection component, the directional odd component, the normal distance perturbation component, and the directional and normal distance interaction component; Trajectory synchronization and antenna coupling correction are performed based on each component, and the results are updated accordingly. A common reflector observation group is established using the updated common reflection component as the stable reflection arrival time; Establish a differential closure relationship for the observation group, determine the pipeline center position, direction and local equivalent radar wave propagation velocity, and obtain candidate propagation velocity values, basic residuals and parameter uncertainties; The observation group is confirmed or split based on the candidate propagation velocity values, each component, and the fundamental residual; Based on each component and the basic residual, supplementary measurement actions are generated, and the survey results are output according to the parameter uncertainty.

2. The method for mapping and identifying underground pipelines using ground-penetrating radar and 3D point cloud technology as described in claim 1, characterized in that, The acquisition of ground-penetrating radar data, 3D point cloud data, attitude data, and odometer data, and the determination of candidate underground reflection areas, includes: Perform time-zero coarse correction on ground-penetrating radar data; The initial reflection trajectory is formed by matching reflection events in adjacent radar channels whose reflection time difference and main phase difference are not greater than preset trajectory time thresholds and trajectory phase thresholds, respectively. The device trajectory is determined based on attitude data and odometer data, and the initial reflection trajectory is projected onto the ground surface represented by the 3D point cloud data according to the initial time index and the device trajectory. The initial reflection trajectories of adjacent survey lines whose projection distance and reflection time difference are not greater than the preset survey line merging distance threshold and survey line merging time threshold, respectively, are merged to form underground reflection candidate areas.

3. The underground pipeline mapping and identification method combining ground-penetrating radar and three-dimensional point cloud as described in claim 2, characterized in that, The step of calculating the local surface normal based on the three-dimensional point cloud data within the candidate region and deploying locally closed measurement units includes: Surface point clouds are extracted from candidate underground reflection regions, and the unit eigenvector corresponding to the smallest eigenvalue of the covariance matrix of the neighborhood point cloud is determined as the local surface normal. The initial extension direction of the candidate underground reflection region is fitted, and the projection of the initial extension direction onto the local tangent plane perpendicular to the local surface normal is determined as the longitudinal unit tangent, and the transverse unit tangent perpendicular to the longitudinal unit tangent is determined. Set up a central measuring point, and set up positive transverse measuring points, negative transverse measuring points, positive longitudinal measuring points, and negative longitudinal measuring points along the transverse unit tangent and the longitudinal unit tangent respectively to form a local closed measuring unit.

4. The underground pipeline mapping and identification method combining ground-penetrating radar and three-dimensional point cloud as described in claim 3, characterized in that, The four-state observation of bidirectional and binormal distances in the measurement unit includes: At the same ground projection location of each measuring point, two different preset normal distances are used to sequentially measure the first normal distance in the forward direction, the first normal distance in the reverse direction, the second normal distance in the forward direction, and the second normal distance in the reverse direction. Record the actual normal distance, attitude data, and odometer data for each state to determine the corresponding device trajectory; The ground-penetrating radar data is resampled according to the three-dimensional point cloud coordinates so that the four states correspond to the same surface projection position.

5. The underground pipeline mapping and identification method combining ground-penetrating radar and three-dimensional point cloud as described in claim 4, characterized in that, The correction of reflection arrival time and orthogonal decomposition yields common reflection components, directional odd components, normal distance perturbation components, and directional and normal distance interaction components, including: Obtain the extrinsic parameters between the pre-calibrated ground-penetrating radar antenna and the attitude sensor; Based on the attitude data, the external parameters, and the local surface normal, determine the phase center of the transmitting antenna, the phase center of the receiving antenna, and the corresponding surface projection point; The air propagation time is calculated based on the distance from each phase center to the corresponding surface projection point and the air wave speed; the reflection arrival time of candidate reflection events under four states in the ground penetrating radar data is extracted as the original reflection arrival time; the air propagation time and the pre-calibrated internal fixed delay of the antenna are subtracted from the original reflection arrival time to obtain the corrected reflection arrival time; Orthogonal decomposition is performed on the arrival times of four corrected reflections at the same location and for the same reflection event to obtain the common reflection component, the directional odd component, the normal distance perturbation component, and the directional and normal distance interaction component.

6. The method for mapping and identifying underground pipelines using ground-penetrating radar and 3D point cloud technology as described in claim 5, characterized in that, The process of synchronizing and correcting the trajectory and antenna coupling based on each component, and then updating it, includes: A synchronization correction model is constructed with the goal of minimizing the root mean square value of the directional odd component and the mileage ratio of reverse measurement, mileage translation, and time synchronization offset as the optimization quantities. The synchronization correction model is solved and spatial resampling and orthogonal decomposition are performed again. Extract direct waves or stable surface reference reflections from a preset time window starting from time zero, determine the time offset corresponding to their normal distance disturbance components as the system group delay, and subtract the system group delay from the arrival time of each corrected reflection. Recalculate and update each component.

7. The method for mapping and identifying underground pipelines using ground-penetrating radar and 3D point cloud combined according to claim 6, characterized in that, The establishment of a common reflector observation group, using the updated common reflection component as the stable reflection arrival time, includes: Retrieve reflection events where the arrival time difference of stable reflections between adjacent measuring points is not greater than a preset continuity threshold; Calculate the principal phase difference between reflection events, extract the effective bandwidth from the spectrum and calculate the relative difference of the effective bandwidth. Establish the correlation between reflection events whose principal phase difference, relative difference of effective bandwidth and interaction components of direction and normal distance are not higher than the associated phase difference threshold, bandwidth difference threshold and interaction threshold of direction and normal distance, respectively. The relative difference of effective bandwidth is the ratio of the absolute value of the difference between two effective bandwidths to their average value. A common reflector observation group is formed by reflection events from at least three non-collinear measuring points; The three threshold values ​​range from 15° to 60°, 0.10 to 0.50, and 0.03ns to 1.00ns, respectively.

8. The method for mapping and identifying underground pipelines using ground-penetrating radar and 3D point cloud combined according to claim 7, characterized in that, The process of establishing a differential closure relation for the observation group, determining the pipeline center location, direction, and local equivalent radar wave propagation velocity, and obtaining candidate propagation velocity values, fundamental residuals, and parameter uncertainties includes: The axis of the candidate pipe segment is represented by the center control point and the directional unit vector of the pipe segment to be estimated. The underground propagation geometric distance is calculated based on the three-dimensional coordinates of each measuring point. The differential propagation time, the horizontal and vertical antisymmetric time difference, and the horizontal and vertical second-order closure quantities are generated by combining the stable reflection arrival time. The candidate values ​​of propagation velocity are obtained from each differential closure relationship. A joint objective function is constructed based on the differential propagation residual, antisymmetric closure residual, second-order closure residual, and the interaction components of direction and normal distance, and weighted robust optimization is performed to determine the pipeline center location, direction, and local equivalent radar wave propagation velocity. The differential propagation residual, antisymmetric closure residual, and second-order closure residual are determined as the basic residuals, and the parameter uncertainty is determined by the parameter covariance matrix.

9. The method for mapping and identifying underground pipelines using ground-penetrating radar and three-dimensional point cloud as described in claim 8, characterized in that, The process of confirming or splitting the observation group based on the candidate propagation velocity value, each component, and the fundamental residual includes: The ratio of the standard deviation to the mean of the candidate propagation speed values ​​is determined as the propagation speed closure index, and the root mean square value of the interaction components of the direction and normal distance within the group is determined as the interaction closure index. When both indicators and each basic residual are not higher than the corresponding threshold, the common reflector observation group is confirmed to be a single linear pipeline. When any index or any basic residual exceeds the corresponding threshold, subgroups are formed according to the principal phase, stable reflection arrival time, spatial projection position, and interaction components of direction and normal distance, and differential closure relationships are established for each subgroup. If no subgroup meets the aforementioned conditions after reaching the preset maximum number of retests, the absolute value of the second-order closure quantity obtained by combining the reflection events of different subgroups is determined as the mixed second-order closure residual. When the mixed second-order closure residual exceeds the diagonal closure threshold, it is marked as a nonlinear underground target group; The propagation speed closure threshold ranges from 0.03 to 0.20, and the threshold for the interactive closure index is the interaction threshold between the direction and the normal distance. The diagonal closure threshold is not less than 0.05 ns.

10. The method for mapping and identifying underground pipelines using ground-penetrating radar and three-dimensional point cloud as described in claim 9, characterized in that, The process of generating supplementary measurement actions based on each component and the basic residual, and outputting surveying results according to parameter uncertainty, includes: The root mean square values ​​of the directional odd component, the normal distance perturbation component, the directional and normal distance interaction component, and the differential propagation residual are successively determined as the directional residual, the normal distance residual, the directional and normal distance interaction residual, and the differential propagation residual index. The larger of the root mean square values ​​of the antisymmetric closed residual and the second-order closed residual is determined as the closed residual. When any residual exceeds the corresponding threshold, the following actions are generated in sequence: reverse retest of the same trajectory, increase normal distance state, increase diagonal measuring point and split common reflector observation group, adjust the center of local closed measurement unit and change measurement baseline. After the supplementary measurements, four-state observations were re-implemented and differential closure relationships were established. The uncertainty of parameters is characterized by the standard deviation of the plane position and the standard deviation of the burial depth. When each closure index, residual and standard deviation does not exceed the corresponding threshold, the axes of adjacent candidate pipe segments whose distance between adjacent endpoints and the angle of orientation are not greater than the preset connection distance threshold and connection angle threshold, respectively, are connected to form the three-dimensional centerline of the underground pipeline. If the aforementioned conditions are still not met after reaching the preset maximum number of retests, an unclosed anomaly flag will be output.

Images