Underground pipeline three-dimensional deformation quantitative inspection method

Abstract

The invention discloses an underground pipeline three-dimensional deformation quantitative inspection method, which comprises the following steps: processing a multi-source sensing sequence to generate a baseline data packet; the method comprises the following steps of: performing first-stage on-line correction, acquiring a medium state agent quantity by applying micro-excitation sampling combined by a micro-amplitude angle and double power to laser, and performing joint inversion to obtain section geometric data subjected to medium correction by utilizing a forward model with endogenous refraction and scattering correction items; second-stage scale robustness processing is carried out, in a unified optimization framework, motor current-based dynamic consistency constraint, geometric conservation constraint and wellhead external reference strong constraint are fused, and a globally consistent pose track is solved; and performing three-dimensional reconstruction and deformation quantification according to the corrected section data and the pose track. According to the invention, the accuracy of geometric measurement of the pipeline and the global consistency of the pose track in a complex environment such as a humid environment can be improved.

Description

Technical Field

[0001] This invention relates to the field of underground pipeline measurement, and specifically to a three-dimensional deformation-based inspection method for underground pipelines. Background Technology

[0002] As the lifeline of urban operations, underground pipelines directly impact the stability and safety of urban water supply, drainage, and gas transmission. With increasing pipeline service life, corrosion, deformation, and damage become more frequent, potentially leading not only to resource leaks and environmental pollution but also major safety accidents such as ground subsidence. Therefore, developing high-precision three-dimensional deformation-based inspection technology for underground pipelines to achieve accurate assessment and early warning of pipeline structural conditions is crucial for ensuring the safety of urban infrastructure, implementing preventative maintenance, and reducing operation and maintenance costs.

[0003] Currently, the internal inspection of underground pipelines has gradually evolved from traditional qualitative, manually interpreted closed-circuit television inspection to quantitative three-dimensional measurement using laser scanning. Existing mainstream technologies typically equip inspection robots with ring laser scanners to acquire the two-dimensional cross-sectional contours of the pipeline frame by frame. This data is then combined with pose information provided by wheel speed odometers and inertial measurement units (IMUs) to stitch together a series of two-dimensional cross-sections into a three-dimensional point cloud model. Some solutions also incorporate Simultaneous Localization and Mapping (SLAM) technology, utilizing visual or geometric features within the pipeline to optimize the pose trajectory, aiming to obtain a more consistent three-dimensional model and achieve quantitative analysis of deformation indicators such as pipeline diameter and ellipticity.

[0004] However, existing technologies face severe challenges in real-world pipeline environments. These challenges are not simply sensor noise issues, but rather stem from deep-seated technical bottlenecks arising from the failure of physical models under complex operating conditions. Current solutions mainly suffer from two problems: first, distortion due to medium disturbances in the measurement reference; and second, instability in the odometry scale of pose estimation. Specifically, existing laser measurement models are generally based on dry, clean air environments. However, in actual humid pipelines, thin water films adhering to the pipe walls, condensed water droplets, or diffuse fog cause refraction and scattering of the laser path, resulting in a systematic, incident angle-related positive offset in ranging. This is not random noise, but a deterministic bias caused by physical model mismatch, which can systematically underestimate the pipe diameter or misjudge an undeformed cross-section as an ellipse. Furthermore, humidity and deposits can cause stick-slip phenomena on robot drive wheels, resulting in a slow and continuous drift in the odometry scale factor based on wheel speed. The scarcity of textures inside the pipe and abrupt changes in illumination further hinder traditional visual SLAM from providing reliable scale constraints. These two problems exacerbate each other, causing the reconstructed 3D model to be distorted locally and inconsistent at the global scale, which severely restricts the reliability of the deformation transformation results. Summary of the Invention

[0005] The purpose of this invention is to provide a three-dimensional deformation-based inspection method for underground pipelines to solve the aforementioned problems existing in the prior art.

[0006] According to one aspect of this application, a three-dimensional deformation-based inspection method for underground pipelines includes:

[0007] Process the raw sequences of multi-source sensors and external mapping benchmark information to generate a baseline sensor data package containing synchronous sensor sequences, wellhead event frame sets, coordinate system and external parameter metadata;

[0008] The baseline sensor data packet is subjected to a first-level online correction to correct for media disturbances, and the geometrically optimized media-corrected cross-sectional data is output.

[0009] By combining the baseline sensor data packet and the medium-corrected cross-sectional data, a second-level scale robustness process is performed to output the global reference pose trajectory.

[0010] Based on the global reference pose trajectory and the cross-sectional data after medium correction, a three-dimensional reconstruction is performed, and a three-dimensional model and deformation quantification index are output.

[0011] Beneficial effects: Through the above technical solutions, the present invention can improve the accuracy of pipeline geometry measurement and the global consistency of position trajectory in complex environments such as humid conditions. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the overall process of a three-dimensional deformable quantitative inspection method for underground pipelines.

[0013] Figure 2 The flowchart illustrates the steps involved in performing a first-stage online correction on the baseline sensor data packet to correct for media disturbances, and outputting geometrically optimized media-corrected cross-sectional data.

[0014] Figure 3 A flowchart of the steps for jointly solving the correction term and cross-sectional geometric parameters.

[0015] Figure 4 This is a flowchart illustrating the steps involved in constructing dynamic consistency constraints.

[0016] Figure 5 This is a flowchart illustrating the logical decision process for determining whether a short time window will ultimately be marked as unrecognizable. Detailed Implementation

[0017] To enable those skilled in the art to better understand the present invention, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, 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.

[0018] The terms "first," "second," etc., used in the specification and accompanying drawings of this invention are used to distinguish different objects, not to describe a specific order. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, apparatus, or product comprising a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or units.

[0019] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of the invention. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0020] Example 1: This example provides a general hardware system and data processing background for the present invention, describing the process of generating baseline sensor data packets. In this example, the three-dimensional deformation-based inspection method for underground pipelines can be executed on an underground pipeline inspection robot system. This robot system integrates one or more ring laser scanners for acquiring the two-dimensional cross-sectional contours of the pipeline; one or more visible light or infrared cameras for capturing the texture and environmental images inside the pipeline; an inertial measurement unit (IMU) containing a three-axis accelerometer and a three-axis gyroscope for measuring the robot's attitude and acceleration; one or more pairs of drive wheels in contact with the inner wall of the pipeline, equipped with wheel speed encoders for measuring the distance traveled; and a current sensor for the drive motor for monitoring the motor's load. All sensor data is timestamped with high precision and is collected and stored by the robot's built-in computing unit.

[0021] According to one aspect of this application, the process of generating a baseline sensor data package involves processing multi-source sensor raw sequences and external mapping benchmark information to generate a baseline sensor data package containing synchronous sensor sequences, a wellhead event frame set, coordinate system, and extrinsic metadata. Specifically, this step, during the data preprocessing and baseline establishment stage before deformation detection and quantization, transforms heterogeneous and asynchronous raw data streams from different sensors into a regularized, aligned base dataset containing necessary metadata, i.e., the baseline sensor data package. This process can be further broken down into the following steps:

[0022] Step one involves performing multi-source time reference alignment. Specifically, this involves reading the original frame sequence of the laser annular cross-section, the original frame sequence of the visible light image, the inertial measurement sequence, the wheel speed and mileage sequence, and the motor current sequence from the stored data stream collected by the inspection robot. Since each sensor has an independent internal clock, its timestamps may have initial deviations and long-term drift. Using a unified system clock reference (e.g., the master clock of the computing unit), and employing timestamp interpolation and drift compensation algorithms, the sampling points of all auxiliary observation sequences (such as IMU attitude, wheel speed, and current) are aligned to the timestamp of each laser annular cross-section synchronization frame, resulting in a frame-by-frame KOM-aligned synchronized sensing sequence. For example, for IMU readings between two laser data frames, linear interpolation or quaternion spherical linear interpolation can be used to estimate its attitude at the laser frame time. This guarantee step ensures that at any given time, there is a state snapshot containing information such as cross-section, attitude, wheel speed, and current, providing a data consistency foundation for subsequent multi-sensor fusion.

[0023] Step two involves screening observation quality and suppressing false points. Specifically, the synchronous sensing sequence may contain abnormal data caused by factors such as strong reflections from inside the tube (e.g., reflections from accumulated water or metal parts), completely black objects (absorbing light and causing no echo), motion blur, or electromagnetic interference. These false points are identified and filtered out by setting a series of thresholds and rules. For example, an upper limit for echo intensity can be set; any laser point exceeding this limit is marked as a strong reflection false point. The edges of shadow areas are identified by detecting spatial discontinuities in the laser point cloud. Motion blur is identified by analyzing drastic morphological changes in the point cloud at the same location between consecutive frames. Identified abnormal data points are removed or marked as low confidence, generating a quality-annotated synchronous frame sequence of the laser annular cross-section, improving the robustness of subsequent processing steps.

[0024] Step 3 involves establishing a coordinate system and solving for extrinsic parameters. Specifically, to unify sensor measurements into a meaningful physical space, multiple coordinate systems need to be established and their transformation relationships determined—the extrinsic parameters. Specifically, this includes defining the equipment's own coordinate system (e.g., with the IMU center as the origin), the laser scanner coordinate system, the camera coordinate system, the pipeline mileage coordinate system (expanded along the pipeline centerline), and an external mapping coordinate system (such as a geodetic coordinate system). Through an offline calibration process, such as simultaneously acquiring laser point clouds and images of a specific calibration board in a calibration field, the relative pose relationships (i.e., extrinsic parameters) between the laser scanner, camera, and IMU can be solved. At the starting point of the inspection task, i.e., the wellhead location, the equipment pose is aligned with known external mapping reference information (such as the GPS coordinates and elevation of the wellhead) manually or automatically, establishing an initial association between the pipeline mileage coordinate system and the external mapping coordinate system. All these coordinate system definitions and extrinsic parameter matrices are stored as coordinate system and extrinsic parameter metadata.

[0025] Step four involves automatically marking short time windows at the wellhead. Specifically, the wellhead is a crucial node connecting the inside and outside of the pipe, and sensor data from its vicinity possesses distinct characteristics that can serve as strong constraints for subsequent global positioning. The original frame sequence of visible light images is analyzed, and a photometric abrupt change detection algorithm (e.g., detecting a sharp increase or decrease in overall image brightness) is used to initially identify the moment the device enters or leaves the wellhead. Furthermore, combined with laser cross-sectional data, the geometric features of the cross-section at this time are evaluated to determine whether they conform to the regular circular or rectangular edges of the well wall. Through this dual verification of photometric and geometric data, a precise time window can be defined, and the data frames within this window are collected into a wellhead event frame set. This set marks the occurrence of wellhead events and is used to trigger subsequent global pose anchor point injection.

[0026] Through the above steps, the raw, messy, multi-source data is processed into a structured baseline sensor data package. This package contains a frame-aligned synchronous sensor sequence, a set of wellhead event frames marked with key locations, and coordinate system and extrinsic metadata defining spatial relationships. This data package forms the basis for all subsequent algorithmic steps.

[0027] Example 2: This example describes the overall process of a three-dimensional deformable quantitative inspection method for underground pipelines according to the present invention, as follows: Figure 1 As shown. The method flow of this embodiment is based on the baseline sensor data packet generated in Embodiment 1, forming a complete technical chain from local correction to global optimization, and then to 3D reconstruction and quantization output. The specific steps are as follows:

[0028] Step S001 involves processing the multi-source sensor raw sequences and external mapping reference information to generate a baseline sensor data package containing synchronous sensor sequences, a wellhead event frame set, coordinate system, and extrinsic parameter metadata. Specifically, in this embodiment, this step, as described in Embodiment 1, provides standardized input data for the entire invention method.

[0029] Step S002 involves performing a first-level online correction on the baseline sensor data packet to correct for media disturbances, outputting geometrically optimized media-corrected cross-sectional data. Specifically, this step addresses the systematic bias in laser ranging caused by media such as thin water films and mist commonly found in underground pipelines. Traditional methods often ignore this influence or use fixed compensation values, leading to distortion in pipeline geometry measurements (such as diameter and ellipticity) in humid areas. This invention uses an active online correction method to treat the media influence as an estimable variable for inversion. Specifically, the system actively applies a small perturbation excitation to the laser scan and establishes a physical forward model that describes the media refraction and scattering effects. Then, an optimization algorithm is used to solve for the media correction parameters and the actual pipeline cross-sectional geometric parameters in real time. The output of this process is a series of media-corrected cross-sectional data whose geometry has minimized the influence of media disturbances, restoring the true contour of the pipeline. Detailed implementation of this step will be discussed in subsequent embodiments.

[0030] Step S003 involves combining the baseline sensor data package and the cross-sectional data after medium correction to perform a second-level scale robustness process, outputting the global reference pose trajectory. Specifically, this step addresses the problem of long-distance pose and odometer drift of the inspection robot caused by sparse features and slippage within the pipeline. Relying solely on wheel speed integration or inertial navigation accumulates non-negligible errors over long distances. This invention constructs a unified optimization framework that integrates various complementary constraint information of different properties to jointly constrain the pose trajectory, achieving long-term scale stability. Specifically, this framework integrates: a global absolute pose anchor point (strong constraint) based on wellhead events and external mapping benchmarks; dynamic consistency constraints based on the motor physical model and differences in wheel speed and inertia; and geometric conservation constraints based on the local geometry of the pipeline maintaining stability in sections with minimal medium influence. By solving this unified optimization problem, a scale-corrected pose trajectory aligned with the global benchmark, i.e., the global reference pose trajectory, can be obtained. Detailed implementation of this step will be discussed in subsequent embodiments.

[0031] Step S004 involves performing 3D reconstruction based on the global reference pose trajectory and the medium-corrected cross-sectional data, outputting a 3D model and deformation quantification indicators. Specifically, using the global reference pose trajectory obtained in the previous step, a series of high-quality, medium-corrected 2D cross-sectional data output in the second step can be accurately placed at their corresponding positions in 3D space. By stitching, meshing, and reconstructing the surfaces of these 2D cross-sections distributed along the trajectory, a detailed 3D model of the pipeline (such as a dense point cloud or triangular mesh model) can be generated. Based on this, the 3D model can be further analyzed to automatically calculate various key deformation quantification indicators, such as the diameter profile varying along the pipeline mileage, the ellipticity profile, and the depth and area of ​​local depressions or protrusions. These indicators provide quantitative and traceable evidence for the structural health assessment of the pipeline.

[0032] As can be seen, the overall method described in this embodiment, through a series of interconnected processes, ensures the accuracy of local cross-section measurements through online correction at the data source and ensures the consistency of global pose trajectory through multi-source constraint fusion. Based on these two aspects, a high-precision 3D model and reliable quantitative indicators are generated, systematically solving the pain points faced by existing technologies.

[0033] Optionally, in an exemplary embodiment, the overall flow of a three-dimensional deformation-based inspection method for underground pipelines may further include: processing the original multi-source sensor sequence to output a synchronous sensor sequence and a set of identified wellhead event frames, and calibrating it in conjunction with external mapping reference information to output a coordinate system and extrinsic metadata; using the synchronous sensor sequence, coordinate system, and extrinsic metadata, performing online correction of medium disturbance to output corrected cross-sectional data; fusing the synchronous sensor sequence, wellhead event frame set, corrected cross-sectional data, and external mapping reference information to perform scale robustness processing to output a mileage scale factor curve and a global reference pose trajectory; and performing three-dimensional reconstruction and deformation based on the corrected cross-sectional data, global reference pose trajectory, and mileage scale factor curve to output a three-dimensional model and deformation indices.

[0034] Understandably, in another embodiment, the overall process of a three-dimensional deformation-based inspection method for underground pipelines can also be as follows: Processing the original multi-source sensor sequence to output a synchronous sensor sequence; identifying the synchronous sensor sequence to output a set of wellhead event frames; calibrating the sequence using external mapping reference information to output a coordinate system and extrinsic metadata; performing online correction of media disturbances on the synchronous sensor sequence, utilizing the coordinate system and extrinsic metadata to output corrected cross-sectional data; performing scale robustness processing on the pose trajectory, fusing the synchronous sensor sequence, the set of wellhead event frames, the corrected cross-sectional data, and the external mapping reference information to output a mileage scale factor curve and a global reference pose trajectory; stitching together the corrected cross-sectional data, utilizing the global reference pose trajectory and the mileage scale factor curve to output a three-dimensional model and deformation-based indices.

[0035] Example 3: This example provides a detailed description of the first-level online calibration steps described in Example 2. This example generates new observables through active micro-excitation sampling, achieving simultaneous inversion and correction of media disturbances (such as refraction caused by thin water films and scattering caused by droplets) and actual pipe geometry within a unified optimization framework. Traditional passive measurement methods cannot decouple the coupling effect between media disturbances and actual geometric deformation, leading to distorted measurement results. This example solves this technical problem through active detection and modeling inversion.

[0036] According to one aspect of this application, the online correction is performed sequentially within a series of short time windows that are consecutive or overlapping along the mileage, wherein each short time window contains a multi-frame synchronous sensing sequence, and the window length is selected to ensure that the pipe geometry and medium state within the window remain quasi-static, so as to meet the data consistency requirement for joint solution.

[0037] In this embodiment, the short time window is a spatiotemporal processing unit. Specifically, it is a data block consisting of a sequence of multiple frames (e.g., 5 to 15 frames) of synchronized laser annular cross-sections continuously acquired along the inspection robot's travel direction (i.e., mileage). The choice of window length is a trade-off: on the one hand, the window needs to be long enough to contain sufficient data frames for effective micro-excitation sampling and statistical inversion; on the other hand, the window must be short enough to ensure that within the few centimeters of pipe section covered by the window, the actual geometric cross-section of the pipe (e.g., diameter, center position) and the physical properties of the environmental medium (e.g., water film thickness, fog concentration) can be considered essentially constant, i.e., satisfying the quasi-static assumption. This assumption is the physical basis for subsequent joint solutions, as it guarantees that all measurement data acquired within the window points to the same actual geometric and medium state to be solved. For example, for a system acquiring data at a frame rate of 10 Hz and traveling at a speed of 0.1 m / s, a short time window containing 10 frames of data corresponds to a mileage of 0.1 meters, which typically satisfies the quasi-static requirement. Optionally, the window length can also be adaptive. For example, in sections where the pipeline geometry or medium condition changes drastically, the system can automatically shorten the window length to improve spatiotemporal resolution.

[0038] Furthermore, a synchronous sensing sequence is extracted from the baseline sensing data packet. This synchronous sensing sequence is then sampled using a micro-excitation method combining micro-amplitude angular oscillation and dual-power laser emission to obtain a medium state surrogate quantity characterizing the medium state. This step is an active information acquisition process designed to address the coupling problem between the medium and geometry. Specifically, the system no longer passively performs measurements in a single mode, but actively applies two small, controlled perturbations—micro-excitations—to the measurement process to elicit different responses from the measured object (pipe and medium).

[0039] Micro-angle oscillation refers to controlling the laser scanner's beam emission direction to make small, periodic angular oscillations in a direction perpendicular to its main scanning plane (i.e., the plane containing the annular cross-section). For example, the oscillation amplitude can be set to ±0.5 degrees. Within a short time window, laser frames can be acquired in a sequence of upward deflection (+0.5 degrees), no deflection (0 degrees), downward deflection (-0.5 degrees), and no deflection (0 degrees). Because the inner wall of the pipe is curved, this tiny emission angle oscillation will cause the laser beam to have an incident angle θ on the pipe wall. inc The angle between the laser beam and the normal to the tube wall surface undergoes a predictable change. This change in the incident angle is crucial for the excitation of the refraction effect, as the ranging error caused by refraction is highly correlated with the incident angle.

[0040] Dual-power laser emission refers to the alternating or paired use of two different laser emission powers within a short time window; for example, a standard power P. stdand a low-power P low (e.g. P) low =0.5*P std Scattering particles (such as fog droplets) in a medium typically respond nonlinearly to laser energy, especially at high powers, which can lead to near-end echo saturation. By comparing echo characteristics at different powers, scattering effects can be identified and quantified more effectively.

[0041] Through the aforementioned combined micro-excitation sampling, the system obtains a series of information-rich raw measurements. Based on this, a medium state surrogate can be constructed, which does not directly measure the physical parameters of the medium itself (such as refractive index), but indirectly characterizes the medium state by calculating some observations sensitive to medium effects. In this embodiment, the medium state surrogate is derived by analyzing the data acquired during micro-excitation sampling, and includes at least one of the following types of physical observations:

[0042] One type involves echo intensity, which characterizes the absorption and scattering properties of a medium, attenuating with a short time window. Specifically, in clean air, the laser echo intensity I roughly follows an inverse square law with distance d. However, the attenuation intensifies in the presence of a scattering medium. Analyzing the change in echo intensity I within a short time window allows for the fitting of a local attenuation coefficient, the magnitude of which is one of the surrogate quantities characterizing the turbidity of the medium.

[0043] Another type is the contour edge energy, which characterizes the edge blurring effect caused by medium scattering. Specifically, ideally, when a laser scans to the edge of a pipe, the distance value changes abruptly. The scattering effect blurs this edge, forming a transition zone. By analyzing the gradient distribution of the contour edge in the range image, its energy or equivalent width can be calculated. The larger this value, the stronger the scattering effect, making it another key proxy for characterizing scattering intensity.

[0044] Based on this, a forward model of the range image with endogenous refraction and scattering correction terms is constructed. Then, using the medium state surrogate quantity, synchronous sensing sequence, coordinate system and extrinsic data, the correction terms and cross-sectional geometric parameters are jointly solved to generate geometrically optimized medium-corrected cross-sectional data.

[0045] This step establishes a mathematical-physical model that connects the actual geometry, medium effects, and the final measurement value; this is the range image forward model. This model accounts for the ranging error e caused by medium perturbations. media Endogenization, as part of the model. An exemplary forward model can be represented as:

[0046] d meas =d geom (P geom )+e media +n=d geom (Pgeom )+f ref (θ inc A ref )+f scat (E blur A scat )+n; where: d meas The distance value actually measured by the laser scanner; d geom The distance from the inspection robot to the inner wall of the pipe is the actual geometric distance, which is determined by the pipe's cross-sectional geometric parameters P. geom (For example, determined by the coordinates of the center of the cross section (cx, cy) and the radius r in the device coordinate system); e media The total ranging error introduced by the medium; n is Gaussian white noise. media It is further broken down into two parts: f ref This is the refraction correction term, which is related to the laser incident angle θ. inc and a set of refraction effect parameters A to be determined ref Related functions; f scat For scattering correction, it is a term related to the contour edge ambiguity E. blur (A medium state proxy) and a set of scattering effect parameters A to be determined. scat Related functions. f ref and f scat The specific function form can be further defined in subsequent embodiments, for example, by using a piecewise linear function or a polynomial function.

[0047] Next, we will perform a joint solution. If we solve the above model directly, since P geom A ref A scat When multiple unknown parameters are coupled together, the problem is ill-conditioned and difficult to obtain a unique solution. Therefore, further, according to one aspect of the invention, such as Figure 3 As shown, the joint solution of the correction term and cross-sectional geometric parameters includes: constructing differential observations to decouple refraction and scattering effects using synchronous sensing sequences acquired at different micro-oscillation angles and different transmission powers. Constructing the differential observations involves: constructing two types of differential observations based on multiple sets of measurement data acquired within a short time window:

[0048] The first category is the angular difference observable Δd. angle Specifically, at the same laser emission power (e.g., both are P), std Under the condition of a slight upward swing angle (corresponding to the incident angle θ), select an angle that swings upward slightly (corresponding to the incident angle θ). inc_up ) and the downward slight swing angle (corresponding to the incident angle θ) inc_down The difference in distance between the two measured laser points is calculated to form an angular difference observable primarily used to separate the effects of refraction. That is: Δdangle_obs =d meas (θ inc_up )-d meas (θ inc_down Since it is within a short time window, the true geometric d corresponding to these two measurement points geom and scattering effect f scat They are considered identical, and therefore cancel each other out after subtraction. Differential observations primarily reflect the change of refraction effect with the incident angle, i.e., Δd. angle_obs ≈f ref (θ inc_up A ref )-f ref (θ inc_down A ref This provides a basis for solving the refraction parameter A. ref It provides strong constraints.

[0049] The second category is power differential observation ΔI. power Specifically, at the same scanning angle (i.e., incident angle θ) inc Under approximately the same conditions, calculate the power P... high With low power P low The difference in echo intensity corresponding to the laser beam forms the power differential observable, primarily used to separate the effects of scattering. That is: ΔI power_obs =I(P) high )-I(P low The difference component affects the scattering parameter A. scat The changes are quite sensitive, in order to solve A scat It provides independent observation information.

[0050] Based on this, a joint optimization objective function is minimized to simultaneously solve for the correction term and the cross-sectional geometric parameters. This joint optimization objective function includes: a fitting residual term matching the differenced observations, and a roundness prior penalty term based on the current cross-sectional geometric parameters. An exemplary joint optimization objective function J can be constructed as: J(P geom A ref A scat )=Σ i w i *(d meas_i -d pred_i ) 2 +w ref *(Δd angle_obs -Δd angle_pred ) 2 +w scat *(ΔI power_obs -ΔI power_pr ed ) 2 +w circ *C(Pgeom ); where: the first term is the fitting residual term of the original distance, d pred_i The first term is the distance prediction calculated based on the forward model; the second term is the fitting residual term from the angle difference observations; the third term is the fitting residual term from the power difference observations; w i ,w ref ,w scat ,w circ These are the weighting coefficients for each item. The fourth item, C(P) geom C(P) is a roundness prior penalty term. Its function is to apply a weak prior constraint that the pipe cross-section is likely to be circular, thereby enhancing the stability of the solution. For example, C(P) geom P can be defined as the sum of squared residuals after fitting all corrected 3D points to a circle. By minimizing the objective function J using nonlinear least squares optimization algorithms such as Levenberg-Marquardt, the optimal P can be solved simultaneously. geom A ref A scat .

[0051] Solving for P geom (Such as the center coordinates and radius) represent the geometrically optimized, medium-corrected cross-sectional data. Compared to the original data, this data more closely approximates the true physical morphology of the pipe, providing high-quality input for subsequent 3D reconstruction and deformation analysis. As an optional implementation, for non-standard circular pipes, the roundness prior penalty term can be replaced with an ellipticity prior or a more generalized smoothness prior. The solution algorithm can also employ more robust M-estimation to suppress the influence of individual outliers.

[0052] Optionally, according to one aspect of this application, a first-level online correction is performed on the baseline sensing data packet to correct for media disturbances, outputting geometrically optimized media-corrected cross-sectional data, such as... Figure 2 As shown, the specific steps can also be as follows: extract the synchronous sensing sequence from the baseline sensing data packet, apply a micro-excitation sampling of the synchronous sensing sequence with a combination of micro-amplitude angle swing and dual-power laser emission, and obtain the medium state surrogate quantity characterizing the medium state; construct a range image forward model with endogenous refraction and scattering correction terms, and use the medium state surrogate quantity, synchronous sensing sequence, coordinate system and extrinsic data to jointly solve the correction terms and cross-sectional geometric parameters to generate geometrically optimized medium-corrected cross-sectional data.

[0053] Example 4: Building upon Example 3, this example describes how to construct an information-driven adaptive closed-loop control system. Example 3 provides a method for media disturbance inversion within a single short time window, but the micro-excitation configurations used (such as angular swing amplitude and power ratio) are preset. However, the media environment within the pipeline (such as water film thickness and mist concentration) is dynamically changing, and a fixed micro-excitation configuration may be suboptimal under certain conditions, or even lead to inversion failure. This example describes an advanced mechanism that can evaluate the inversion quality online and dynamically optimize subsequent micro-excitation configurations accordingly, elevating the measurement process from open-loop execution to closed-loop self-optimization.

[0054] According to one aspect of this application, a first-level online correction is performed on the baseline sensing data packet to correct for medium disturbances, and geometrically optimized medium-corrected cross-sectional data is output. The method further includes: evaluating the correction parameter uncertainty of the correction term obtained by joint solution and the cross-sectional geometric parameters; and adaptively updating the micro-excitation configuration file for subsequent sampling periods based on the evaluated correction parameter uncertainty, wherein the micro-excitation configuration file controls the angular swing amplitude, power ratio or sampling window length of the micro-excitation sampling.

[0055] Specifically, this step is performed after each joint optimization solution in Example 3 is completed. The uncertainty of the corrected parameters is evaluated. Mathematically, the uncertainty of a parameter is inversely proportional to the sharpness of its objective function J at the optimal solution. A sharp valley (high curvature) indicates that the parameters are well constrained by the data, resulting in low uncertainty; conversely, a flat valley (low curvature) indicates that the data constrains the parameters weakly, resulting in high uncertainty. A preferred method for quantifying uncertainty is to utilize the Hessian matrix H, a byproduct of the optimization process, i.e., the objective function J with respect to all parameters (P... geom A ref A scat The second-order partial derivative matrix of is given. Near the optimum, the covariance matrix Cov of the parameters can be approximated by the inverse of the Hessian matrix, i.e., Cov≈H. -1 Each correction parameter (e.g., refractive parameter A) ref A certain component a in ref_k The uncertainty σ(a) ref_k σ(a) can be estimated by taking the square root of the diagonal elements of the covariance matrix. ref_k = sqrt(Cov(k,k)) where: σ(a ref_k ) is parameter a ref_k The standard deviation of the covariance matrix is ​​its uncertainty; Cov(k,k) is the k-th diagonal element of the covariance matrix Cov.

[0056] Furthermore, based on the assessed uncertainty of the correction parameters, an adaptive update is performed. The system compares the calculated uncertainty σ with a preset confidence threshold T. conf The comparison is performed, and the update logic is executed based on the comparison result.

[0057] When the uncertainty is high (e.g., σ(a) ref_k )>T conf_ref This indicates that the data generated by the current micro-excitation configuration is insufficient to identify the refraction effect. At this point, the system will update the micro-excitation configuration profile to strengthen the corresponding excitation in the next short time window. Specifically, the angular oscillation amplitude can be increased, for example, from ±0.5 degrees to ±0.8 degrees. A larger oscillation amplitude will produce a larger range of incident angle variation, in the angle difference observation Δd angle_obs This generates a stronger signal, which helps to more accurately identify the refractive parameter A. ref This reduces the uncertainty in the next solution.

[0058] When the uncertainty of the scattering parameters is high, similarly, if the scattering parameter A... scat If the uncertainty is too high, the system can adjust the power ratio, for example, by increasing the high and low power P. high and P low The difference between them, to enhance power differential observation ΔI power_obs The signal-to-noise ratio.

[0059] When the overall uncertainty is high or the solution fails to converge, it may indicate that the signal-to-noise ratio of the current environment is extremely low. In this case, the system can selectively increase the sampling window length, for example, increasing the number of laser frames included in the window from 10 to 15. A longer window provides more measurement data points, which helps to suppress noise and improve the stability of the solution through statistical averaging.

[0060] In this way, the system forms a complete feedback control loop: measurement-inversion-uncertainty assessment-adjustment of measurement strategy-new measurement. This loop enables the inspection method to intelligently adapt to various complex and unknown pipeline environments.

[0061] Furthermore, the micro-amplitude angle swing used in the micro-excitation sampling has an angle swing control sequence that is pre-generated based on identifiability sensitivity analysis. This generation process includes: predicting the identifiability sensitivity of different scanning incident angles to refraction and scattering correction terms by performing small perturbation simulation of the measurement process; and selecting the angle in the high sensitivity range based on the sensitivity prediction result to construct the angle swing control sequence.

[0062] In addition to the aforementioned real-time adaptive control, this embodiment also includes a method for prior optimization design of the micro-excitation strategy to ensure its basic performance. Specifically, for the design of the angle swing sequence, a identifiability sensitivity analysis can be performed in advance. This analysis aims to identify which incident angles θ... inc The combination of these parameters is useful for identifying refractive parameter A. ref Most sensitive or efficient. The specific steps are as follows:

[0063] Step 1: Perform small-perturbation simulation, specifically, based on the forward model d established in Example 3. pred =d geom +f ref (θ inc A ref )+…+n, set a nominal pipe geometry and medium parameters in the simulation environment.

[0064] Step two, calculate the sensitivity. Specifically, the sensitivity is defined as the observed quantity (here, the angle difference observation Δd). angle_pred ) relative to the parameter to be determined (A) ref The rate of change of θ, i.e., the partial derivative S(θ) inc_up ,θ inc_down )=|&(Δd angle_pred ) / &A ref |(& is for partial derivative). The larger this value, the smaller the parameter A is. ref The larger the change in parameter A, the greater the change in the observed value. ref The easier it is to be identified.

[0065] Step three involves constructing an optimal sequence. Specifically, by calculating the sensitivity across all possible combinations of incident angles, a sensitivity map can be obtained. The system will then select the angle combinations located in the high-sensitivity region to construct the actual angle swing control sequence, ensuring that the selected angles are within the mechanical motion range of the robot hardware. For example, analysis might indicate that the sensitivity is highest when one incident angle is around 5 degrees and another around 15 degrees; in this case, the angle swing amplitude should be designed to keep the incident angles within this range.

[0066] Based on this, the dual-power laser emission used in micro-excitation sampling has its high and low power pairing dynamically determined according to the echo saturation state. This determination process includes: analyzing the historical echo intensity distribution in the recently acquired synchronous sensing sequence to predict the echo saturation threshold online; and generating a pair of high and low power emission based on the predicted saturation threshold. This pairing is to ensure that the high-power echo does not overflow the dynamic range of the sensor, and the low-power echo has a distinguishable signal-to-noise ratio.

[0067] This step involves further dynamic optimization of the dual-power excitation strategy. The varying materials, colors, and moisture levels of the pipe's inner wall result in significant variations in its laser reflectivity. A fixed high- and low-power pairing may work well in one section, but in another, it may lead to high-power echo saturation (signal overflow) or an excessively low low-power echo signal-to-noise ratio (signal submerged in noise). Specifically:

[0068] The system maintains a sliding histogram of echo intensity I containing data from the most recent hundreds of frames.

[0069] The laser sensor's ADC (analog-to-digital converter) has a maximum reading, such as 4095. If the histogram shows a large number of echo intensity values ​​concentrated at this maximum value, it indicates echo saturation. The system can then dynamically assess the current risk of saturation based on this.

[0070] Based on saturation risk assessment, the system dynamically adjusts the high- and low-power transmission pairing for the next short time window. high ,P low Its adjustment logic is as follows:

[0071] If signs of saturation are detected, reduce P appropriately. high The value is such that its echo intensity falls just near the upper limit of the dynamic range, for example, the target intensity is 95% of the maximum value.

[0072] P low The value needs to be related to P high Maintain a suitable ratio (e.g., P) low =0.5*P high This ensures that the echo intensity is much higher than the noise floor of the sensor, thus guaranteeing a identifiable signal-to-noise ratio.

[0073] As an alternative, the adaptive update logic described above can be implemented without using a simple threshold rule, but instead employ a more complex control algorithm, such as a PID controller or a fuzzy logic controller, to achieve smoother and more robust policy adjustments.

[0074] This embodiment transforms the micro-excitation sampling process from a fixed open-loop operation into an intelligent, closed-loop, and self-optimizing cognitive sensing process through real-time uncertainty feedback, prior sensitivity analysis, and dynamic power spectrum adjustment.

[0075] Example 5: This example elevates and unifies the media correction parameters obtained from discrete window-by-window solutions in Examples 3 and 4 into a single field quantity continuously distributed along the pipeline mileage and circumference, i.e., the equivalent media thickness distribution. This distribution map is then used to provide feedback guidance throughout the measurement and estimation process. The media distribution within the pipeline (e.g., water accumulation is mainly at the bottom) exhibits strong spatial continuity and correlation. Integrating this prior information into the estimation process yields a smoother, more robust, and physically more clearly defined media distribution map, providing crucial weighting information for subsequent processes. The specific steps are as follows:

[0076] Step S401 involves constructing an equivalent medium thickness distribution that continuously varies along the pipeline mileage and circumferentially by introducing a joint optimization process that incorporates spatial continuity and geometric consistency priors for multiple medium-corrected cross-sectional data generated along the pipeline mileage.

[0077] Specifically, this step will take the series of discrete correction parameters A output from Examples 3 and 4. ref and A scat These parameters are obtained as input observations. Although they are solved independently within each short time window, they collectively reflect a latent, continuously varying physical field. This step describes the inversion of this latent physical field.

[0078] Specifically, a two-dimensional mesh is defined to parameterize this equivalent medium thickness distribution field. The two dimensions of this grid are the pipeline's mileage (m) and the cross-sectional circumferential angle. The value T of each node on the grid ij This represents the distance in m. i Circumferential angle The equivalent medium thickness or degree of influence at the location.

[0079] Next, construct a grid with all grid node values ​​{T} ij Let} be the joint optimization objective function for the variables. This function typically includes the following terms:

[0080] Data fitting term J data This requirement is to reverse the scene. This best explains the discrete correction parameters previously solved within each short time window. In other words, it... After substituting into the physical model, its predicted correction effect should be consistent with A. ref and A scat The observed values ​​match.

[0081] Spatial continuity prior (smooth prior) J smooth This term introduces a universal law of the physical world: the distribution of a medium usually changes smoothly, without drastic, discontinuous jumps. This term is typically modeled as a field. The sum of squares of gradients in the mileage and circumferential dimensions. For example, J smooth =μ*Σ((T) i +1,jT i ,j) 2 +(T i ,j+1-T i ,j) 2 Here, μ is the weight that controls the smoothing intensity. This term acts like a low-pass filter, effectively suppressing noise caused by single-window solution errors.

[0082] Geometric consistency prior J geom This term introduces more advanced prior knowledge. For example, in a straight pipe section identified as continuous and undamaged, its true geometric curvature should be constant. If the medium correction is imperfect, it can cause undue fluctuations in the calculated geometric curvature. This term penalizes such geometric inconsistencies caused by medium residuals.

[0083] By minimizing the overall objective function J field =J data +J smooth +J geom This allows for the calculation of the values ​​of all grid nodes, resulting in a complete and continuously varying equivalent medium thickness distribution map. This map visually demonstrates which areas within the pipe (such as the bottom) are severely affected by the medium. This distribution map is not only the final analysis result but also crucial intermediate information, which can be used for macroscopic feedback adjustments to the entire methodology, improving overall performance.

[0084] Step S402, based on the equivalent medium thickness distribution, performs at least one of the following feedback adjustment operations: one operation is that feedback modulation is used for the micro-excitation configuration profile in subsequent sampling cycles to enhance sampling in areas where the medium's influence is significant. This is a region-based adaptive sampling strategy. When the inspection robot is about to enter an area predicted to be wet based on the generated distribution map (i.e.,... When the value is high in a certain area, the system can proactively and in advance switch to the micro-excitation configuration profile of the enhanced sampling mode. For example, it can automatically increase the amplitude of the angle swing and the power ratio, and even increase the density of sampling points. Compared with the window-by-window passive feedback in Embodiment 4, this feedforward adjustment method has a faster response and can better cope with regional changes in the medium environment.

[0085] Another operation is to dynamically adjust the weight of geometric observations during scale robustness processing to reduce the contribution of data from the medium-affected region to the scale solution. This is another manifestation of the high coupling of information across multiple modules in this invention. In the subsequent scale robustness processing (detailed in later embodiments), an important source of constraint is the geometric conservation constraint, which assumes that the diameter of the pipe is constant in certain sections. However, in the wet region where the medium influence is significant, despite the correction in Embodiment 3, the accuracy and reliability of its cross-sectional geometric measurements may still be lower than in the dry region. In this case, the equivalent medium thickness distribution map becomes an excellent confidence map. When constructing the unified optimization objective function for scale robustness, for each geometric conservation constraint, the weight of the observations is adjusted according to the location m. A value is used to assign a dynamic weight. Specifically, The higher the value of a region, the lower its weight for geometric observations. This dynamic weighting mechanism effectively prevents measurement errors remaining in wet areas from contaminating the global-scale solution, thus improving the robustness and accuracy of the final pose trajectory.

[0086] As an alternative, the solution for the equivalent medium thickness distribution can be obtained without global one-time optimization. Instead, a sequential estimation method based on Kalman filtering or particle filtering can be used to update the field distribution segment by segment to meet the needs of online real-time processing.

[0087] Example 6: This example provides an alternative or supplementary implementation method to enhance real-time performance. While the joint optimization inversion method described in Example 3 boasts high accuracy and robustness, its iterative solution nature places certain demands on computational resources. In some high-speed inspection or embedded computing scenarios with limited resources, it may be difficult to guarantee frame-by-frame real-time processing. Therefore, this example describes a fast media correction method based on offline learning and online retrieval fine-tuning. By trading space for time, it reduces the complexity of online computation while maintaining correction accuracy.

[0088] According to one aspect of this application, in order to improve real-time performance, the process of generating the media-corrected cross-sectional data may further include an enhancement step, specifically, using a media state proxy quantity to retrieve initial local mapping prototype parameters from a preset local reversible mapping library that stores the mapping relationship between media state and correction parameters; and based on the local mapping prototype parameters, performing online linearization fine-tuning on the current cross-sectional data to quickly generate an enhanced version of the media-corrected cross-sectional data.

[0089] Specifically, the complex mapping relationship between medium states and optimal correction parameters, while difficult to express analytically with simple formulas, exhibits strong local similarity. In other words, similar medium states (characterized by medium state surrogate variables) typically correspond to similar optimal correction parameters. This step utilizes this characteristic to transform the complex nonlinear optimization problem into a fast table lookup and linear interpolation problem.

[0090] For example, this step requires building a locally reversible mapping library offline. This library is a pre-learned and populated structured database, such as a KD-tree, hash table, or trained neural network model. The library stores a large number of (state-parameter) pairs.

[0091] The state vector S, specifically, serves as a query key and consists of a series of surrogate variables calculated within the current short time window that characterize the state of the medium. For example, S = [p decay ,p blur ,p intensity_mean ,...], where p decay It is the echo intensity attenuation coefficient, p blur It refers to the blurriness of the outline edge, p intensity_mean It refers to the average echo intensity, etc.

[0092] The parameter vector P, as the value returned by the query, specifically contains the optimal correction parameters solved using the complete joint optimization method of Example 3 in this state. For example, P = [A ref A scat The library can be constructed by calibration in a laboratory or by distilling and filling it with the fine-mode algorithm of Example 3 on a large amount of real pipeline inspection data. Each processed short time window contributes one (S,P) data point. During online inspection, when media calibration is required for a new short time window, the following operations are performed:

[0093] State vector calculation and retrieval, specifically, according to the method in Example 3, calculate the medium state vector S of the current window. current With S current To query the index, a search is performed in the locally reversible mapping library. The goal of the search is to find the data stored in the library that corresponds to S. current The closest one or more state vectors S lib Proximity here can be defined using Euclidean distance or other metrics. The retrieval returns S. lib The corresponding parameter vector P lib The P lib These are the initial local mapping prototype parameters. This can be understood as P... lib It is the optimal correction parameter among historical experiences most similar to the current situation.

[0094] Furthermore, online linearization fine-tuning is performed; specifically, due to S current With S lib There will always be a slight difference between them, ΔS = S current -S lib Use P directly lib This will introduce a certain approximation error. To compensate for this error, a quick fine-tuning is needed. Assume that in S... lib Within the neighborhood of P, the mapping relationship between parameter P and state S is approximately linear, i.e., ΔP ≈ J * ΔS. Here, J is the Jacobian matrix of parameter P with respect to state S, J = &P / &S (& represents the partial derivative), which characterizes the sensitivity of the correction parameter to state changes. This Jacobian matrix J can also be calculated offline and stored in the library along with (S,P). Therefore, the online fine-tuning process simplifies to a simple linear update operation: P final =P lib +J lib *(S current -S lib Wherein: P final These are the correction parameters used after fine-tuning; P lib and J lib The prototype parameters and corresponding Jacobian matrices are retrieved from the library. This calculation process involves only one matrix-vector multiplication and one vector addition, and its computational cost is much less than the iterative optimization in Example 3.

[0095] Furthermore, the obtained correction parameter P final Substituting the range image forward model described in Example 3, for all the original range values ​​d within the current window... meas After calibration, enhanced media-calibrated cross-sectional data can be quickly generated.

[0096] As an optional implementation, the fast correction method of this embodiment can be used in combination with the fine correction method of Embodiment 3. For example, the system defaults to using the fast method of this embodiment for real-time processing, while running the full optimization of Embodiment 3 in the background at a lower frequency. When the correction residual predicted by the fast method is large, or the state vector S current With S in the library lib When the distance exceeds a threshold (indicating the encounter of unknown new operating conditions), the system can seamlessly switch to fine mode to ensure accuracy in the worst-case scenario. The calculation results in fine mode can be used to dynamically and online update and expand the local invertible mapping library, enabling the system to continuously learn and evolve.

[0097] Example 7: This example provides a detailed description of the scale robustness process. This step aims to address the inherent odometer drift and scale uncertainty issues in long-distance pipeline inspections. Examples 3 to 6 ensured the accuracy of local two-dimensional cross-sectional measurements, while this example stitches these accurate local measurements together in a globally consistent three-dimensional space aligned with the external world. To this end, this example constructs a unified optimization framework that integrates the advantages of multiple information sources (including external benchmarks, dynamic models, and geometric constraints).

[0098] Specifically, the second-level scale robustness process is achieved by constructing and solving a unified optimization framework that integrates at least two types of constraints: dynamic consistency constraints based on the driving torque derived from the motor current sequence and the attachment prior; and geometric conservation constraints extracted based on the geometric stability of the local cross section, wherein the weight of the constraint is dynamically modulated by the equivalent medium thickness distribution.

[0099] In this embodiment, scale robustness is used to solve for the precise pose (position and attitude) of the inspection robot at every moment, specifically to correct the scale factor drift accumulated in its odometer integral. This scale factor *s* defines the proportional relationship between the odometer increment calculated by the wheel speedometer or inertial navigation and the actual odometer increment in the physical world; ideally, it is 1, but it can change due to factors such as slippage and wheel wear. This invention uses a unified optimization framework to measure the complete pose trajectory T of the inspection robot. traj (Including position, attitude, and scale factor s(t) that change over time) as variables to be optimized, and introducing multiple constraints to jointly determine its optimal solution.

[0100] Furthermore, the unified optimization framework also utilizes strong external benchmark constraints based on wellhead events. The generation and injection of these strong constraints include: automatically identifying and selecting key wellhead frames from the set of wellhead event frames based on the consistency assessment of photometric mutations and wellbore geometric edges; and generating pose anchors with global coordinates by binding the poses of the key wellhead frames with external mapping benchmark information. These pose anchors are then injected into the unified optimization framework as hard constraint terms.

[0101] This step provides a global reference point for the entire relative drift pose trajectory. The wellhead is the only location where in-tube measurements can be directly correlated with ground mapping coordinates (such as GPS or total station points).

[0102] Within the set of wellhead event frames generated in Example 1, the system further searches for the frame or frames with the richest and most reliable information as wellhead keyframes. The selection is based on the consistency assessment of photometric abrupt changes and wellbore geometric edges. Specifically, the system searches for a frame whose visible light image brightness gradient reaches a local peak (corresponding to the moment of most drastic change in light intensity at the wellhead), and whose simultaneous laser cross-section data can be matched with a standard wellbore geometric model (such as a circle or rectangle) with extremely low residual. This dual condition ensures that the selected keyframe is not only temporally accurate but also geometrically reliable.

[0103] For a selected wellhead keyframe, the system binds its pose to pre-stored external mapping reference information (such as the precise geodetic coordinates of the wellhead) and solves for it. This is a small-scale nonlinear optimization problem to solve for the precise pose (position and orientation) of the keyframe in the global coordinate system. This solved pose with global coordinates constitutes the pose anchor point. In the unified optimization framework, this pose anchor point is introduced as a hard constraint or a strongly constrained term with extremely high weight. This means that the optimized global pose trajectory T traj At the moment it passes the wellhead keyframe, its pose must be exactly equal to or very close to the value of that pose anchor point. This effectively eliminates global position and heading drift throughout the trajectory.

[0104] In addition to the aforementioned strong constraints, the framework also incorporates at least two types of soft constraints that run throughout the entire trajectory:

[0105] One is the dynamic consistency constraint, which specifically utilizes the robot's inherent physical motion laws. It is independent of external environmental features, making it particularly important in sections with sparse texture inside the pipe. The specific construction of this constraint will be detailed in subsequent embodiments.

[0106] Another constraint is the geometric conservation constraint, which specifically utilizes the geometric properties of the pipe itself. It assumes that in an ideal pipe section without deformation or media interference, the pipe diameter measured from high-quality cross-sectional data should remain constant. The construction of this constraint and how its weights are modulated by the equivalent media thickness distribution generated in Example 5 will also be detailed in subsequent examples.

[0107] As an alternative implementation, the unified optimization framework can be constructed as a factor graph model. In this graph, the robot pose to be optimized (including the scale factor) is the variable node, while the wellhead anchor constraints, dynamic gas consistency constraints, geometric conservation constraints, etc., mentioned above serve as factor nodes connecting the corresponding variable nodes. The entire scale robustness problem is then transformed into solving the maximum a posteriori probability estimation problem of this factor graph, which can be solved efficiently using mature optimization libraries such as g2o and CeresSolver.

[0108] Example 8: This example provides a detailed description of dynamic consistency constraints. This example provides a continuous, intrinsic constraint source for pose estimation in long-distance, featureless environments, independent of external visual or geometric features. It is rooted in fundamental physical principles of robotics: in the absence of slippage, the distance traveled by integrating wheel speed should be consistent with the displacement integrated by the inertial measurement unit; and the occurrence of slippage is closely related to the driving torque applied by the robot and the current surface adhesion coefficient. This example translates this physical insight into mathematical constraints.

[0109] According to one aspect of this application, such as Figure 4 As shown, the construction of dynamic consistency constraints includes: quantifying the differential residual between wheel speed integral mileage and inertial integral mileage; mapping the motor current sequence to drive torque via the motor physical model; and establishing an objective term to penalize the differential residual, wherein the weight of the penalty is dynamically adjusted by the mapped drive torque.

[0110] Specifically, the construction of this constraint consists of the following three parts:

[0111] The first part is the differential residual between the quantized wheel speed integral mileage and the inertial integral mileage. The motion of the inspection robot can be measured by at least two independent internal sensors, specifically:

[0112] Wheel speed integral mileage d wheel (t) is obtained by time integration of the wheel angular velocity measured by the wheel speed encoder. This measurement is directly related to the scaling factor s(t) to be optimized, d wheel (t)=∫s(t)*v wheel (t)dt. Where v wheel (t) is the wheel speed encoder reading. d wheel The advantage of s(t) is that it has high accuracy and no zero bias at low speeds, but its fatal weakness is that it cannot detect slippage, and its scale factor s(t) will drift slowly due to wheel wear, air pressure changes, or even slight slippage.

[0113] Inertial integral mileage d inertial (t) is obtained by performing two time integrations on the acceleration measured by the IMU (after gravity compensation and attitude rotation to the world coordinate system). That is, d inertial (t)=∫∫R(t)*(a raw (t)-g)dtdt. Where a raw Here, R is the raw acceleration reading, R is the attitude matrix, and g is the gravity vector. The advantage of inertial measurement is that it directly measures the actual motion of the body, unaffected by wheel slippage. However, its disadvantages include zero bias and noise, and quadratic drift can occur during long-term integration.

[0114] Difference residual rodom (t) represents the difference between these two at every time t: r odom (t)=d inertial (t)-d wheel (t). This residual r odom The physical root of the problem is either that the scale factor s(t) is inaccurate or that actual physical slippage has occurred. This constraint is designed to intelligently distinguish between these two possibilities.

[0115] The second part maps the motor current sequence to driving torque using a physical model of the motor. Motor current is an excellent and easily measurable indicator of the drive system's state. For the DC motors commonly used in inspection robots, its output torque τ... drive (t) and the current I flowing through the armature motor There exists an approximately linear relationship between τ and t: drive (t)=K t *I motor (t)-τ friction Wherein: K t This is the motor's torque constant, a parameter that can be found in the motor manual or obtained through simple calibration; I motor (t) is the current value measured in real time by the current sensor on the motor drive board; τ friction It is the equivalent frictional torque within the system, which can be calibrated as a constant or a function of velocity. This model allows for real-time, non-intrusive knowledge of how the robot moves.

[0116] The third part is to establish the objective term for penalizing the difference residuals. This objective term J dyn This will be added to the unified optimization objective function in Example 7, and its basic form is J. dyn =Σw dyn (t)*(r odom (t)) 2 It is important to note the penalty weight w. dyn (t) is not a fixed constant, but is affected by the driving torque τ drive The dynamic regulation of (t). This is based on a profound physical insight, namely the slippage sparsity prior, which is described in the following steps.

[0117] In one exemplary embodiment, the objective term for penalizing the differential residuals is configured as a function to apply a slippage sparsity prior; the function applies a strong penalty to non-zero differential residuals to force them to converge to zero when the drive torque is low, and relaxes the penalty when the drive torque is high to allow physically reasonable slippage to occur.

[0118] This step ensures that this constraint can intelligently identify the source of error. The aforementioned slippage sparsity prior means that, in most cases, the robot moves stably, and slippage is a low-probability, sparsely occurring event. The physical prerequisite for slippage is that the driving torque approaches or exceeds the maximum static friction between the wheels and the road surface. Based on this, the weighting function w dyn (t) is designed to be a function of the driving torque τ drive The nonlinear function of (t), specifically:

[0119] When the driving torque is low, the robot is simply moving smoothly and at a low speed, far from reaching its adhesion limit, and physically, slippage is impossible. Therefore, the observed non-zero residual r odom All residuals (t) should be attributed to the inaccurate scaling factor s(t). To force the optimizer to adjust s(t) to eliminate this residual, the weights w... dyn (t) is set to a large value, and a strong penalty is applied.

[0120] When the driving torque is high, the robot may be climbing a slope, accelerating, or traversing a very slippery surface. The driving torque is approaching the adhesion limit, and slippage is physically possible. The residual r at this point... odom The residual s(t) likely reflects the actual physical slippage, and the scale factor s(t) should not be forced to absorb this residual. Therefore, the weight w dyn (t) is set to a smaller value, which relaxes the penalty and allows the optimizer to believe that this is a genuine slippage event.

[0121] An example weight function w dyn (t) can be constructed in the following form: w dyn (t)=W max / (1+exp(k*(τ drive (t)-τ limit )));wherein: W max It is a very large constant weight; τ limit This represents the estimated driving torque under the current road surface adhesion limit; k is the coefficient of the control function steepness. It can be seen that this is an inverse S-shaped function, implementing the aforementioned control logic.

[0122] In an exemplary embodiment, the attachment prior in the dynamic consistency constraint is generated by: mapping the driving torque calculated from the motor current sequence into a preset road condition experience library that characterizes various road surface physical characteristics; and dynamically assigning an attachment prior value to the current working condition based on the mapping result.

[0123] This step resolves the key parameter τ in the aforementioned weighting function. limitThe problem involves the dynamic determination of the adhesion limit τ. The road surface conditions within the pipeline are constantly changing (dry concrete, slippery PVC, cast iron pipes with deposits, etc.), thus the adhesion limit τ... limit It is not static either. Specifically:

[0124] An offline database is built to store the physical properties of different road surface types. For example: {Type: "Dry Concrete", Coefficient of Adhesion μ: 0.8}, {Type: "Wet and Slippery PVC", Coefficient of Adhesion μ: 0.2}, {Type: "Muddy Deposits", Coefficient of Adhesion μ: 0.3}.

[0125] When running online, the system can determine the current road surface type in several ways. One method is to use visible light camera images for texture analysis and material classification. Another, more direct method is to analyze the behavior of the drive system itself. For example, the system can continuously monitor the drive torque τ. drive (t) and actual acceleration a inertial The relationship between (t) is established. If a difference between wheel speed and inertial speed is observed under a very small driving torque (i.e., slippage occurs), the system can determine that it is currently on a road surface with a low coefficient of adhesion. By matching these online observation features with an experience base, the system can dynamically assign the most probable prior value of adhesion coefficient μ to the current operating condition. prior .

[0126] μ was determined prior The adhesion limit torque τ can then be calculated. limit =μ prior *F normal *R wheel F normal It is the supporting force on the wheels (related to the robot's weight and the pipe's inclination angle), R wheel It's the wheel radius. This dynamically updated τ limit It was then used in the weighting function.

[0127] As an alternative, the motor model can be further considered to account for back electromotive force and viscous friction to obtain a more accurate torque estimate. Weighting function w dyn (t) can also be expressed in other robust statistical functions, such as the Huber function or Tukey's biweight function, to further enhance robustness to abnormally large slip events.

[0128] This embodiment combines the physical model of the motor, kinematic differences, and dynamic attachment priors to construct a strong, adaptive dynamic constraint, providing support for achieving scale robustness of pose trajectories under various complex working conditions.

[0129] Example 9: This example provides a detailed description of geometric conservation constraints. If the dynamic constraints of Example 8 are based on the robot's introspection, then the geometric constraints of this example are based on the understanding of the external environment, i.e., the geometry of the pipe itself. As a man-made structure, the pipe's cross-section (e.g., diameter) should remain constant in sections where no deformation occurs. Utilizing this prior knowledge, another strong constraint source that positively complements the dynamic constraints can be provided for the scale factor s(t) of the pose trajectory. This example does not blindly trust all geometric observations, but rather uses the equivalent medium thickness distribution map generated in Example 5 to dynamically and intelligently weight the reliability of the geometric observations.

[0130] According to one aspect of this application, the construction of geometric conservation constraints and their weight modulation include: identifying geometrically stable sections from the medium-corrected cross-sectional data along the mileage to form candidates for geometric conservation constraints; querying the equivalent medium thickness distribution; and assigning a dynamic weight to each candidate geometric conservation constraint based on the degree of medium influence in the section where the candidate geometric conservation constraint is located, wherein the weight is negatively correlated with the degree of medium influence.

[0131] Specifically, the construction and application of this constraint involves the following two key steps:

[0132] The first step involves identifying geometrically stable sections from the medium-corrected cross-sectional data along the mileage, forming candidates for geometric conservation constraints. This step automatically identifies the pipe sections most likely to satisfy the geometric conservation premise from the continuous cross-sectional data stream.

[0133] Specifically, the input data is a series of cross-sectional geometric parameters generated in Example 3 or 6 after media correction, such as the cross-sectional diameter D(m) varying with mileage m.

[0134] The system employs a sliding window approach, traversing the diameter sequence D(m) along a mileage m. Within each window, the statistical variance Var(D) of the diameter is calculated. If this variance is less than a preset stability threshold T... stable If the geometry of the pipe segment within this window is considered stable, it can be considered a candidate for geometric conservation constraints. For example, in a straight pipe segment composed of standard pipe sections, the diameter variance will be very small; while in the deformation zone or at the interface, the diameter variance will increase.

[0135] For each identified stable segment (e.g., from mileage m) start to m end The corresponding geometric conservation constraint can be formalized into an objective term J. geom The objective is to penalize the diameter D(m) of each point within the segment compared to the median diameter D within the segment. median The deviation. That is: J geom =Σm_end m_start w geom (m)*(D(m)-D median ) 2 It is important to note that the measured diameter D(m) is directly related to the scale factor s(m) in the pose trajectory. If the scale factor s(m) drifts (e.g., increases), all measured distances will increase accordingly, causing the calculated diameter D(m) to systematically increase as well. Therefore, by minimizing J... geom The optimizer is driven to adjust the curve of the scale factor s(m) so that the corrected diameter sequence remains as flat as possible within the stable region, thus achieving scale constraint.

[0136] The second step involves querying the equivalent medium thickness distribution and assigning a dynamic weight to each candidate geometric conservation constraint based on the degree of medium influence within the segment containing the constraint. This step demonstrates the information flow and coupling between different modules of the invention. Not all geometrically stable segments are equally reliable. The reliability of a stable diameter sequence measured in a severely flooded segment is far lower than that measured in a dry segment, because the former may be affected by incompletely corrected medium residuals.

[0137] For each of the above candidate geometric conservation constraints, the system will query the equivalent medium thickness distribution map generated in Example 5. By considering the mileage range [m] of this constraint start ,m end By integrating or averaging the influence of the medium within a given area (T value), the overall humidity score for that section can be obtained. wet .

[0138] Target item J geom The weight w in geom (m) is no longer a fixed constant, but is set to be related to Score wet Negative correlation. That is, the higher the degree of influence of the medium (Score) wet The larger the value of the geometric constraint, the greater its influence (weight w) within the unified optimization framework. geom The lower the weight, the better. An example weight function could be: w geom =W base *exp(-γ*Score wet ); where: W base γ is the base weight; γ is the coefficient that adjusts the composite sensitivity of the weights. This dynamic weighting mechanism ensures that the scaling process can intelligently trust geometric observations from wet and slippery sections while giving greater trust to geometric observations from dry and clean sections. This avoids the risk of global scaling being contaminated by local medium residuals, and improves the final accuracy and reliability of scaling robustness.

[0139] As an alternative, the identification of stable sections can be based not only on the variance of the diameter, but also on higher-order geometric features such as the ellipticity and curvature of the cross-section. In addition to querying the degree of influence of the medium, the weighted modulation function can also integrate other quality assessment indicators, such as laser echo intensity and inversion uncertainty (from Example 4), to form a more comprehensive confidence assessment system.

[0140] Example 10 describes a proactive observability enhancement step within the method system of this invention to address extremely challenging scenarios. While the unified optimization framework described in Examples 7-9 provides robust pose estimation in most cases, its final accuracy still depends on the continued effectiveness of the input constraints. In certain highly challenging sections, such as a long, geometrically perfect, featureless pipe covered with a uniform and stable water film, resonant states may occur where dynamic constraints are weakened due to smooth motion, geometric constraints are weakened due to the lack of morphological features, and media correction residuals are difficult to discern, leading to the accumulation of pose estimation uncertainty within that section. This example describes a closed-loop control method that can automatically diagnose and proactively resolve such weak observability problems.

[0141] According to one aspect of this application, a three-dimensional deformation-based inspection method for underground pipelines further includes a trajectory reinforcement step triggered on demand, comprising: outputting a report evaluating the quality and uncertainty of the global reference pose trajectory, and automatically identifying weakly observable sections based on the report; when a weakly observable section is identified, triggering the device to perform a short-distance reverse movement covering the section to generate bidirectional consistent observation constraints; adding the newly generated bidirectional consistent observation constraints to a unified optimization framework for scale robustness processing for secondary optimization, and outputting an updated global reference pose trajectory.

[0142] Specifically, this on-demand trajectory reinforcement step transforms the entire inspection process from a linear data collection-processing-output flow into an intelligent process with self-assessment and proactive intervention capabilities.

[0143] Specifically, the system outputs a report evaluating the quality and uncertainty of the global baseline pose trajectory, and automatically identifies weakly observable segments based on this report. During the solution process of the unified optimization framework in Example 7 (e.g., factor graph-based optimization), not only can the optimal estimate of the pose trajectory be obtained, but its information matrix or covariance matrix can also be analyzed to obtain the uncertainty of the pose estimate at each time step.

[0144] Specifically, generating a quality and uncertainty report involves the system calculating the covariance matrix Cov(m) of the pose for each mileage point m in real-time or near real-time. This report can be understood as an uncertainty curve that varies along the mileage; for example, the uncertainty U(m) of the overall pose at that mileage point can be quantified using the trace or the largest eigenvalue of the covariance matrix.

[0145] The automatic identification of weakly observable sections specifically involves the system setting an uncertainty performance threshold T. uncertainty and a segment length threshold L weak When the system detects that the uncertainty U(m) is consistently higher than T... uncertainty The mileage length exceeds L weak When (for example, when the uncertainty is high for 5 consecutive meters), this section [m] start ,m end This is automatically identified and marked as a weakly observable segment. Physically, this means that within this segment, all available constraints (dynamics, geometry, etc.) are insufficient to constrain the robot's pose drift within acceptable limits.

[0146] Furthermore, when a weakly observable section is identified, the device is triggered to perform a short-distance reverse movement covering that section, generating bidirectional consistent observation constraints. Once the triggering condition is met, the system's decision module sends a high-level instruction to the robot's motion control module.

[0147] The robot will smoothly stop moving forward and then begin a controlled, short-distance reverse journey along the original path, with the goal of completely and again scanning the newly marked weakly observable segment [m]. start ,m end During reverse travel, all sensors (laser, IMU, wheel speedometer, etc.) functioned normally and recorded data.

[0148] Now, for the same physical pipe segment [m] start ,m end The system possesses two independent measurement datasets: one from the forward movement and the other from the reverse movement. These two datasets, operating within the same global coordinate system, should describe the same static 3D scene. The bidirectional consistent observation constraint is the mathematical expression of this physical fact. Specifically, the system extracts the local point cloud PC generated by the forward scan. fwd and the local point cloud generated by reverse scanning PC bwd The optimal relative pose transformation ΔT between the two point clouds is calculated using a point cloud registration algorithm (e.g., the Iterative Nearest Point (ICP) algorithm or a variant thereof). fwd_bwd This ΔT fwd_bwdIt is a newly generated, strong geometric constraint that directly relates the pose on the forward trajectory segment to the pose on the reverse trajectory segment.

[0149] Based on this, the newly generated bidirectional consistent observation constraints are added to the unified optimization framework of scale robustness processing for secondary optimization, and the updated global reference pose trajectory is output.

[0150] This newly calculated relative pose transformation ΔT fwd_bwd This factor is added as a new factor (or constraint) to the factor graph model in Example 7. This factor connects one or more pose variable nodes on the forward trajectory with the corresponding pose variable nodes on the reverse trajectory, and its function is similar to a local, high-precision loop closure detection.

[0151] After the new constraint is added, the system will trigger a global optimization solution again. Since this new constraint provides strong geometric closed-loop information for the originally sparse and weakly observable segment, the optimizer can use it to significantly eliminate the accumulated uncertainty in the segment and accurately correct the two trajectories (forward and backward) of the segment as well as the adjacent trajectories.

[0152] After the second optimization, the output is the updated global reference pose trajectory. Compared to the initially generated trajectory, this updated trajectory significantly improves the accuracy and reliability of the previously identified weakly observable segments.

[0153] As an alternative, the decision to trigger trajectory reinforcement can not only rely on pose uncertainty but also combine other quality indicators, such as the residual magnitude of constraint terms and the convergence speed of the optimization process, to form a multi-dimensional decision system. Besides reverse movement, for robots with lateral movement capabilities, observability can also be actively stimulated by performing small-range lateral translations. This embodiment endows the inspection system with the ability to reflect and remedy, transforming the traditional open-loop data acquisition mode into a closed-loop active perception mode that guarantees the quality of the final result and executes on demand.

[0154] Example 11: This example specifically describes one of the final output stages of the method flow of the present invention: how to attach a traceable, spatialized quality and uncertainty report to the generated 3D model and quantitative indicators. Traditional pipeline inspection deliverables are often bare 3D models or deformation data tables, leaving users unable to know which areas of the model have high accuracy or which areas' measurements may be significantly affected by environmental factors. This example solves this problem by systematically collecting, transmitting, and mapping quality-related metadata throughout the entire processing flow, achieving transparency and traceability of the final results.

[0155] According to one aspect of this application, a three-dimensional deformation-based inspection method for underground pipelines further includes the step of generating a traceable quality and uncertainty report for the three-dimensional model. This step includes: extracting the equivalent medium thickness distribution and wellhead constraint quality records generated during the processing as sources of uncertainty for assessing the final result; mapping and labeling these two sources of uncertainty onto corresponding spatial segments of the three-dimensional model to form a segmented confidence assessment; and combining the confidence assessment with deformation-based indicators to generate a quality and uncertainty report.

[0156] Specifically, the process of generating this report involves projecting quality control information from each stage of the method of this invention onto a three-dimensional model.

[0157] Specifically, the equivalent medium thickness distribution and wellhead constraint quality records generated during the processing are extracted as sources of uncertainty in the final results. This invention inherently generates various intermediate results characterizing data quality within the processing flow, eliminating the need for additional quality checks. This step only requires extracting them as the raw input for uncertainty assessment.

[0158] Equivalent medium thickness distribution maps are the best source for assessing uncertainties introduced by medium disturbances. (See figure) The higher the value, the more severe the influence of the medium on the original laser measurement in that area. Even after calibration, the likelihood and magnitude of residual error are still higher than in clean, dry areas. Therefore, It can be directly used as a medium-dependent uncertainty field.

[0159] Wellhead constraint quality records are the source for evaluating the spatial distribution of global positioning accuracy. In regions near high-quality wellhead anchors (with small solution residuals and good geometric matching), the accuracy of the global pose is strongly constrained and therefore very high. As the mileage moves away from this anchor, the pose uncertainty gradually increases due to drift accumulation until the next high-quality anchor is encountered. Therefore, the quality of each wellhead constraint and its position on the mileage axis together define the spatial distribution of global positioning accuracy.

[0160] Furthermore, these two sources of uncertainty are mapped and labeled onto corresponding spatial segments of the 3D model, forming segmented confidence assessments. This step transforms the aforementioned abstract quality information into intuitive and visual annotations on the 3D model.

[0161] For any triangular facet or point on a 3D model, it has unique spatial coordinates, from which its mileage m and circumferential angle can be retrieved.

[0162] Furthermore, for each point v in the model, a comprehensive confidence score Conf(v) can be calculated. This score can be a weighted fusion of multiple factors. An exemplary calculation method is as follows:

[0163] Where: Conf(v) is the overall confidence level of point v (the higher the value, the more reliable it is); Here, is the mileage and circumferential angle corresponding to point v; f(T) is a function negatively correlated with the equivalent medium thickness T, the larger T is, the lower the score for this item; g(...) is a function negatively correlated with the distance Dist to the nearest high-quality anchor point and the mass Q of that anchor point. anchor Positively correlated functions; w media and w loc These are the weights of two sources of uncertainty.

[0164] The calculated confidence score Conf(v) can be stored as an additional attribute of the vertices of the 3D model (e.g., vertex color). In visualization software, the confidence distribution of the entire pipeline model can be visually displayed using pseudo-color rendering. For example, green represents high-confidence areas (dry and clean, close to high-quality wellheads), yellow represents medium-confidence areas, and red represents low-confidence areas (severely slippery, far from all wellheads).

[0165] Based on this, the confidence level assessment is combined with deformation quantification indicators to generate a quality and uncertainty report. The final delivered report is not just a list of deformations, but a comprehensive report that combines what it is and how credible it is.

[0166] Each deformability indicator listed in the report (e.g., a depression with a diameter of -5% detected at kilometer marker 105.3) is accompanied by a confidence score or interval derived from the segmented confidence assessment of its location. For example: ...depression, the overall confidence score for this area is 95% (high); or, ...ellipticity is 1.1, with a 95% confidence interval of [1.08, 1.12].

[0167] The report will also include links to the original evidence that led to the change in confidence level. For example, for a deformation alarm with low confidence, users can click to view the equivalent medium thickness distribution map of that section, or the uncertainty curve of the global trajectory, to understand the specific reason for the low confidence level (e.g., the alarm is located in a severely flooded area).

[0168] As an alternative implementation, the uncertainty fusion model can employ more complex probabilistic models, such as Gaussian processes, to more accurately model the spatial propagation of uncertainty. The final report can also be an interactive 3D scene where users can freely navigate and select any location to query detailed quantitative indicators and corresponding multi-source quality information. This embodiment enhances the engineering value and decision support capabilities of the invention's technological achievements by establishing a complete quality traceability chain from processing to final results.

[0169] Example 12: This example provides a concrete, reproducible numerical calculation case for the joint inversion method of medium perturbation based on micro-excitation sampling. This example quantifies how the present invention utilizes active micro-excitation sampling, particularly micro-amplitude angular oscillations, to generate differential observations that can be used for decoupling, and accurately inverts the true geometric parameters from the original measurement data contaminated by medium perturbation.

[0170] Suppose an inspection robot is inspecting a horizontal pipe with a nominal inner diameter of 800 mm. What is its actual geometric inner radius r? true The distance is constant at 400.0 mm. The rotation center of the robot's ring laser scanner coincides with the geometric center of the pipe (i.e., there is no eccentricity). There is a uniform air density gradient inside the pipe caused by temperature difference. This gradient produces a refraction effect on the laser path similar to a thin lens, resulting in an overestimation of the distance measured.

[0171] To simplify the explanation, this case study will focus on analyzing the correction process for the refraction effect, and will make the following assumptions:

[0172] Geometric parameter to be determined: estimated radius r of the pipe est .

[0173] The medium parameter to be determined: A, a single equivalent parameter describing the refraction effect. ref The unit is mm / rad (millimeters per radian).

[0174] Forward model simplification: A simplified linear refraction model is adopted, i.e., the ranging error e ref The incident angle θ of the laser beam on the tube wall inc The angle between the laser beam and the normal to the tube wall is proportional to: e ref =A ref *θ inc The final measured value d meas For d meas =r est +A ref *θ inc +n, where n is the measurement noise.

[0175] Micro-excitation configuration: The system performs two consecutive measurements aligned with the side of the pipe (e.g., at the 3 o'clock position), using micro-amplitude angular oscillations upwards α. up = +1.0° and downward α down = -1.0° (relative to the horizontal plane).

[0176] Specifically, a ground truth value is set to simulate the physical world, and the goal of the algorithm is to solve for this value.

[0177] Defined true geometric radius: r true =400.0mm.

[0178] Set actual medium parameters: A ref_true = 20.0 mm / rad.

[0179] Next, the true incident angle of the laser is calculated in these two micro-excitation measurements. For a scanner located at the center with radius r, a laser beam with an angle α to the horizontal plane is emitted into a circular pipe of radius r. Based on geometric relationships and the sine theorem, its incident angle θ on the pipe wall can be accurately calculated. inc It equals α. Therefore:

[0180] The angle of incidence θ when swinging upwards inc_up =1.0°.

[0181] The angle of incidence θ when swinging downwards inc_down = 1.0°. (Note: The angle of incidence is the angle with the normal. In this symmetrical geometry, the angles of incidence are the same whether the angle is tilted upwards or downwards by the same amount, but for generality, they are still distinguished.) To more generally represent the change in angle, assume that due to slight misalignment, the actual angle of incidence is θ. inc_up =1.0° =0.01745 rad; θ inc_down =1.1° =0.01920 rad.

[0182] The true ranging error caused by the medium can be calculated:

[0183] Error e of upward measurement ref_up =A reftrue *θ inc_up =20.0 * 0.01745 = 0.349 mm.

[0184] Error e of downward measurement ref_down =A reftrue *θ inc_down =20.0 * 0.01920 = 0.384 mm.

[0185] Generate a simulated, contaminated raw measurement distance and add a small amount of measurement noise (e.g., ±0.05 mm):

[0186] The distance d measured upwards meas_up =r true +e ref_up +n1=400.0+0.349+0.03=400.379mm.

[0187] The distance d measured downwards meas_down =r true +e ref_down +n2=400.0+0.384-0.04=400.344mm.

[0188] In practical applications, d meas_up and d meas_down The direct output of the inspection system sensors is the input data that the algorithm of this invention deals with. It can be seen that if these contaminated data (for example, taking their average value (400.36 mm)) are used directly to estimate the radius, a system error of about 0.36 mm will be generated.

[0189] Based on the above information

[0190] d meas_up =400.379mm,d meas_down = 400.344mm, θ inc_up =0.01745rad,θ inc_down =0.01920 rad, and r est and A ref_est It is unknown.

[0191] According to Example 3, the difference in measured distances is calculated:

[0192] Δd angle_obs =d meas_up -d meas_down =400.379-400.344=0.035mm.

[0193] According to the forward model, theoretically Δd angle_pred =(r est +A ref_est *θ inc_up )-(r est +A ref_est *θ inc_down ). Due to r est The cancellation, Δd angle_pred =A ref_est *(θ inc_up -θ inc_down Let the observed value be approximately equal to the predicted value Δd. angle_obs ≈Δd angle_pred A can be solved directly. ref_est :

[0194] A ref_est =Δd angle_obs / (θ inc_up -θ inc_down = 0.035 / (0.01745-0.01920) = 0.035 / (-0.00175) = -20.0 mm / rad. (The negative sign here indicates the relationship between the defined angle and the error increment, with the same physical meaning). It can be seen that the solved medium parameter -20.0 mm / rad is almost identical in value to the set true value of 20.0 mm / rad, proving the ability of the differential method to identify medium parameters.

[0195] Based on the above steps, a quantitative model of the medium effect (i.e., A) was obtained. ref_est =-20.0), can be used to calculate and eliminate errors in the original measurement:

[0196] Estimated upward measurement error e ref_est_up =A ref_est *θ inc_up =-20.0*0.01745=-0.349mm.

[0197] Estimated downward measurement error e ref_est_down =A ref_est *θ inc_down =-20.0*0.01920=-0.384mm.

[0198] Furthermore, the original measurements are corrected:

[0199] Corrected upward distance d corr_up =d meas_up -e ref_est_up =400.379-(-0.349)=400.03mm (d should be the diameter here) corr_up =d meas_up -e ref_est_up =400.379-0.349=400.03mm). (This should be d) corr_up =d meas_up -e ref_up_est =400.379-(20.0*0.01745)=400.379-0.349=400.03mm, based on the solved A ref (Correct the absolute value). Recalculate:

[0200] A ref_est =Δd angle_obs / (θ inc_up -θ inc_down = 0.035 / (0.01745-0.01920) = -20.0;

[0201] e ref_est_up =-20.0*0.01745=-0.349mm;

[0202] d corr_up =d meas_up -e ref_est_up =400.379-(-0.349)=400.728mm. The calculation is incorrect; the logic needs to be re-examined.

[0203] Forward model: d meas =r true +A ref *θ inc ;

[0204] Measured value: d meas_up =r true +A ref *θ inc_up d meas_down =r true +A ref *θ inc_down ;

[0205] Difference: d meas_up -d meas_down =A ref *(θ inc_up -θ inc_down Solve for A. ref :

[0206] A ref_est =(d meas_up -d meas_down ) / (θ inc_up -θ inc_down = (400.379 - 400.344) / (0.01745 - 0.01920) = 0.035

[0207] -0.00175 = -20.0 mm / rad. The solution is correct. Correction:

[0208] r est_from_up =d meas_up -A ref_est *θ inc_up =400.379-(-20.0*0.01745)=400.379+0.349=400.728mm. Still incorrect.

[0209] The problem is that, physically, A ref It should be a positive value (a positive angle of incidence results in a positive error). However, the above calculation yielded a negative value because θ... inc_up <θ inc_down And d meas_up>d meas_down This does not match the model. Next, we will correct the simulation data to make it conform to the model. Assume θ inc_up =1.1° =0.01920 rad; θ inc_down =1.0° =0.01745 rad;

[0210] e ref_up =20.0 * 0.01920 = 0.384 mm; e ref_down =20.0 * 0.01745 = 0.349 mm;

[0211] d meas_up =400.0+0.384+0.03=400.414mm; d meas_down =400.0+0.349-0.04=400.309mm.

[0212] Recalculate, input: d meas_up =400.414mm; d meas_down = 400.309 mm; θ inc_up =0.01920 rad; θ inc_down = 0.01745 rad;

[0213] Solve A ref A ref_est =(400.414-400.309) / (0.01920-0.01745)=0.105 / 0.00175=60.0mm / rad. This deviates from the true value of 20.0 due to noise. We will temporarily ignore the noise and use the noise-free measurement value: d meas_up_noiseless =400.384mm,d meas_down_noiseless =400.349mm;

[0214] A ref_est = (400.384 - 400.349) / (0.01920 - 0.01745) = 0.035 / 0.00175 = 20.0 mm / rad. The solution is consistent with the actual value.

[0215] Perform geometric correction.

[0216] r est_from_up =d meas_up_noiseless -A ref_est *θ inc_up =400.384-20.0*0.01920=400.384-0.384=400.0mm.

[0217] r est from down =dmeas down noiseless -A ref est *θ inc down =400.349-20.0*0.01745=400.349-0.349=400.0mm.

[0218] By correcting the two measurements separately, two consistent height estimates of the true radius were obtained, both being 400.0 mm. The average of these estimates can be taken as the final geometric parameter estimate, r. est =400.0mm.

[0219] By comparison, it can be found that the traditional method (direct averaging of measured values) does not estimate the medium parameter (i.e., treats it as 0), and the estimated pipe radius is 400.3665 mm (obtained by averaging two measurements of 400.384 mm and 400.349 mm), which has an error of +0.3665 mm compared with the true value of 400.0 mm; the medium parameter A estimated in the method of this invention... ref_est The value is 20.0 mm / rad, the estimated pipe radius is 400.0 mm, and the error between the actual value of 400.0 mm and the actual value is 0.0 mm.

[0220] The method of this invention improves the accuracy of pipe radius estimation by estimating the medium parameter (20.0 mm / rad), achieving zero error; while the traditional method, which does not consider the medium parameter, results in a certain error in the estimation result.

[0221] This numerical calculation example demonstrates that traditional measurement methods introduce systematic geometric measurement errors in the presence of medium disturbances. However, this invention, by actively applying a micro-amplitude angular oscillation and constructing an angular difference observable, can successfully measure the unknown medium influence parameter A. ref The geometric parameters *r* are separated and precisely solved. Based on this, the original measurement data is compensated and corrected to obtain geometric measurement results consistent with the physical true values. This fully verifies the effectiveness and advancement of the method of this invention in improving the accuracy of three-dimensional quantitative inspection of underground pipelines. In real-world scenarios containing noise, by jointly optimizing a large number of points within a short time window, the influence of noise can be effectively smoothed, obtaining equally accurate estimation results.

[0222] Example 13: This example details a series of key strategy initializations and model building steps performed before executing the media correction algorithm. These details are crucial for ensuring the stability, robustness, and correctness of the physical meaning of the correction algorithm.

[0223] Specifically, a sampling strategy initialization file is generated. At the start of an inspection task or during a system cold start, a reasonable initial operating point, i.e., an initial micro-excitation configuration file, needs to be provided for the adaptive micro-excitation control system. The generation of this file is based on the inherent physical properties and task constraints of the inspection robot platform, as detailed below:

[0224] Read the technical manual or calibration data of the laser scanner servo mechanism to determine its maximum permissible angular velocity, angular acceleration, and physical limits of mechanical rotation stroke. Based on this, set a safety upper limit for angular sway with a safety margin. For example, if the physical limit is ±3°, the initial sway amplitude can be conservatively set to ±0.5° to ensure that it will not damage the equipment.

[0225] Evaluate the battery capacity of the inspection robot and the rated power consumption of the laser. High-power P in a dual-power laser emission strategy. high When choosing a device, it is necessary to consider not only the signal-to-noise ratio, but also to ensure that its average power consumption during long-term operation is within the overall energy budget, so as to avoid the mission's battery life falling short of expectations due to local high-intensity sampling.

[0226] Based on the robot's typical inspection speed and the expected frequency of pipeline geometric changes, an initial short time window length is set. For example, for conventional municipal pipelines, the initial window can be set to include 10 frames of laser data, which is usually sufficient to achieve a good balance between satisfying quasi-static assumptions and ensuring the amount of solution data. The initial angle swing amplitude, initial power ratio, initial window length, and other parameters obtained from the above evaluation and setting are collectively stored as an initial micro-excitation configuration file for subsequent adaptive control modules to load and use at the start of the task.

[0227] Furthermore, the key independent variable (incident angle) is precisely calculated in the forward model d of Example 3. meas =d geom +f ref (θ inc A ref In the equation () + ... + n, the laser incident angle θ inc It is a key, known independent variable. The accuracy of its calculation directly affects the accuracy of the medium parameter inversion. The calculation process of this incident angle is as follows: For any measurement point p obtained by laser scanning... scan (In the scanner coordinate system) using the extrinsic parameters calibrated in Example 1, the coordinates are transformed to the device body coordinate system. Based on a preliminary pipe cross-section fitting result within the current short time window (e.g., performing a fast least-squares circle fitting on the original point cloud), the approximate center C of the cross-section can be obtained. approx and radius r approx . Through measurement point p scan and the center C approxA radius line can be determined, and the direction of this radius line is point p. scan The approximate normal direction of the pipe wall n wall The direction of the laser beam, v beam It is known (determined by both the scanner's scanning angle and its slight oscillation angle). Incident angle θ inc That is, these two vectors n wall and v beam The angle between them can be calculated using the dot product formula: θ inc =acos((v beam ·n wall ) / (|v beam |*|n wall This calculation process will accurately calculate the incident angle for all laser points participating in the joint optimization within each short time window.

[0228] Based on this, the form of the constraint correction term function is specified. To enhance the stability and convergence of the joint solution in Example 3, it is necessary to modify the refraction correction term f in the forward model. ref and scattering correction term f scat The specific function form should be reasonably set and constrained.

[0229] Refraction correction term f ref (θ inc A ref Specifically, based on physical laws and experimental observations, the refractive effect of media such as thin water films typically increases monotonically with increasing incident angle. Therefore, a preferred functional form is to set it as a function of the incident angle θ. inc The function is a monotonically piecewise linear function. For example, its nonlinear response can be approximated by a three-segment linear function with two inflection points. For the parameter set A describing this function... ref (i.e., the slope and inflection point of each segment) impose feasible region constraints. For example, based on historical statistics and physical consistency, the slope must be positive and less than a certain upper limit to avoid physically meaningless solutions during the optimization process.

[0230] Scattering correction term f scat (E blur A scat Specifically, similarly, the scattering effect typically varies with the turbidity of the medium (as determined by the proxy E). blur The gain increases with the increase of E, but saturation also exists. Therefore, its functional form can be set as a monotonic function with saturation suppression. For example, a sigmoid function or a gain-compression function can be used, which has a gain of saturation suppression. blur It grows approximately linearly at lower levels, and at E blur At very high values, it gradually approaches a plateau to simulate the saturation characteristics of the scattering effect. For the parameter set A describing this function...scat Similar feasible region constraints are also imposed.

[0231] The refined definition and constraints on the sampling strategy, key independent variables, and model function form provided by the above-mentioned approach lay the foundation for the successful implementation of the correction algorithm. This is a key step in the present invention's transition from a theoretical model to a robust engineering implementation.

[0232] Example 14: This example provides crucial, quantifiable input information for adaptive control and quality reporting, and specifically describes online quality monitoring and evaluation during media correction. A robust industrial-grade algorithm not only needs to process data but also needs introspection capabilities—that is, the ability to know the quality of its current solution in real time. This example describes the mechanism that endows the invention with this introspection capability, enabling the system to quantify its own uncertainties and identify conditions where the model may fail.

[0233] Specifically, quantifying the parameter uncertainty, this step, after the joint optimization process in Example 3 converges, yields a set of optimal geometric and dielectric correction parameters P. opt =[P geom A ref A scat Then execute. Assign a confidence interval or error bar to each component of this set of optimal solutions.

[0234] As described in Example 4, the uncertainty of parameter estimation and the optimization objective function J at the optimal point P opt The valley's sharpness (i.e., curvature) is inversely proportional to the valley's depth. Directly calculating this curvature (i.e., the Hessian matrix H) can be computationally expensive. In nonlinear least squares problems, there exists an efficient approximation method that utilizes the Jacobian matrix J already calculated during the optimization process. res To approximate the information matrix F (the information matrix can be considered as an approximation of the Hessian matrix), we use the Jacobian matrix J. res The information matrix F is defined as the matrix of first-order partial derivatives of all residual terms (including distance fitting residuals, difference observation residuals, etc.) with respect to all parameters to be determined. F ≈ J res T *W*J res ; where: J res T W is the transpose of the Jacobian matrix; W is a diagonal weight matrix containing the weights w of each residual term of the objective function in Example 3. i ,w ref This information matrix F quantitatively describes the total amount of information contributed by all measurement data to constrain each parameter.

[0235] According to statistical theory, the covariance matrix Cov(P) of the parameter to be determined can be obtained by inverting the information matrix F: Cov(P) ≈ F -1 The covariance matrix Cov(P) is a square matrix, and the k-th element on its diagonal is Cov(P). [k,k] That is, the kth parameter P k Variance Var(P) k Therefore, the uncertainty of this parameter (usually represented by the standard deviation σ) can be calculated: σ(P) k = sqrt(Var(P) k ))=sqrt(Cov(P) [k,k] Through this process, the system can calculate each parameter, such as the pipe radius r and the refraction parameter a. ref_k Specific, quantified uncertainties were calculated, for example, σ(r) = 0.15 mm. This assessment of the uncertainty of the correction parameters forms the data basis for all subsequent quality-related decisions.

[0236] In addition to uncertainty assessment based on the optimization results, the system also needs to directly identify some adverse operating conditions that may lead to model failure from the raw data level. This step uses joint decision rules to mark those unidentifiable windows where the inversion results are unreliable.

[0237] Plateau detection (i.e., echo saturation detection) is used to identify situations where the laser sensor's ADC (analog-to-digital converter) reaches its maximum range. In practice, the system analyzes the echo intensity values ​​of all laser points within a short time window. If it finds that more than a certain percentage (e.g., 5%) of the points have echo intensity values ​​precisely equal to the sensor's maximum reading (e.g., 4095 for a 12-bit ADC), the system raises a saturation flag. This is because once saturation occurs, the true information about the echo intensity is lost, which severely compromises the effectiveness of scattering correction models that rely on intensity information and power differential observations.

[0238] Whether a short time window is ultimately marked as unidentifiable or low-confidence depends on a multi-condition logical judgment rule. A preferred rule is, for example... Figure 5 As shown, IF (saturation flag is TRUE) OR (uncertainty of all key parameters σ(P)) k Exceeding its corresponding performance threshold T k THEN marks the current window as unrecognizable. ENDIF, where the performance threshold T... k It is set based on experience or task requirements, for example, T radius It can be set to 1mm, T ref_param It can be set to 50%.

[0239] Once this unidentifiable window marker is generated, it will have two important downstream effects: First, it will act as one of the strongest signals to trigger the adaptive control module in Embodiment 4, forcing it to adjust the sampling strategy in the next window (e.g., if it is marked due to saturation, the laser power must be reduced); Second, the marker will be passed down as metadata along with the cross-sectional data of the current window, and the 3D model of that section will be explicitly marked as low confidence in the quality report of Embodiment 11.

[0240] As an alternative, uncertainty assessment can also employ Monte Carlo methods or bootstrap methods, estimating the parameter distribution by resampling the residuals. This approach is more robust to non-Gaussian noisy models but incurs greater computational cost. The decision rule for the unidentifiable window can also be a more complex decision tree or a classifier trained by a machine learning model, incorporating more quality-related features.

[0241] Example 15: This example specifically describes the quality control of key constraints in the method of the present invention, as well as the refined processing and quantification of the final three-dimensional delivery results.

[0242] Key constraints (wellheads) form the basis of global localization in Example 7, and their quality determines the global accuracy of the 3D model. Therefore, rigorous quality assessment of each potential wellhead constraint is crucial. After the system initially determines a wellhead keyframe using photometric and geometric features and solves its global pose, it will jointly evaluate it from the following dimensions and calculate a comprehensive quality score Q. anchor (between 0 and 1):

[0243] The first dimension is geometric matching confidence, which evaluates the degree to which the laser cross-section point cloud of the keyframe matches the standard wellbore geometry (circular or rectangular). A high-quality match should have an extremely low root mean square (RMS) fitting residual. The magnitude of this residual can be normalized and used as a score.

[0244] The second dimension is the solution residuals and convergence, which evaluates the final residual size and the number of iterations required for convergence of the nonlinear optimization solver when bound to an external mapping benchmark. Low residuals and fast convergence generally imply good constraint consistency.

[0245] The third dimension is the observation geometry, which assesses the observation geometry conditions during measurement. For example, grazing observations of the wellbore due to severe robot tilt should be avoided, as such observations have very weak constraints on depth and normals. The angle between the line-of-sight angle and the wellbore normal should be close to 0 degrees.

[0246] Taking all the above factors into account, a weighted quality score Q can be obtained. anchor This wellhead constraint quality record will be stored and used in subsequent unified optimization.anchor The higher the anchor value, the greater the weight of the corresponding constraint term. It is also a key input for evaluating the spatial distribution of global positioning uncertainty in Example 11.

[0247] After obtaining the scale-robust global pose trajectory and the high-quality cross-sectional point cloud sequence after medium correction, these discrete point clouds need to be transformed into a continuous, complete, watertight three-dimensional surface model (i.e., a mesh model).

[0248] Specifically, odometer coordinate registration and point cloud stitching are performed. Using the final global pose trajectory, the cross-sectional point cloud of each frame is transformed to the same global coordinate system. Through accumulation and stitching, a complete and preliminary dense point cloud is formed.

[0249] Further, surface reconstruction and mesh generation are performed. A surface reconstruction algorithm is used to recover the underlying surface topology from the disordered dense point cloud. Preferred algorithms include Poisson surface reconstruction, which has good smoothing effects on point clouds with some noise; or the moving cube algorithm, which is highly efficient. The output of this step is an initial 3D triangular mesh model.

[0250] Based on this, hole repair and topology optimization are performed. Due to sensor occlusion, false point removal, or poor local data quality, the initial mesh may contain some small holes or topological errors. The system automatically detects these holes and uses boundary flow-based interpolation algorithms or Laplace smoothing methods to reliably repair small holes, generating a topologically complete and watertight mesh. This step is crucial for subsequent accurate volumetric or hydraulic analyses.

[0251] After obtaining a refined 3D mesh model, the system will automatically and quantitatively extract the deformation quantification indicators that users are most concerned about.

[0252] One key metric is the diameter versus ellipticity profile. The extraction process involves virtually slicing a 3D mesh along the mileage axis at fixed, high-density intervals (e.g., every 5 cm), resulting in a series of 2D cross-sectional profiles. For each profile, a robust elliptic fitting algorithm is used to calculate its major axis, minor axis, and azimuth. This yields the diameter versus mileage profile (typically represented by the minor axis or equivalent diameter) and the ellipticity versus mileage profile (the ratio of the major and minor axes) along the entire pipeline mileage. These are crucial metrics for assessing radial deformation (such as flattening) of the pipeline.

[0253] Another metric is the estimation of second-order curvature and local depressions / protrusions. To detect more refined local defects (such as pits and obstacles), the system calculates the second-order curvature (e.g., Gaussian curvature and mean curvature) of each vertex on the 3D mesh surface. Pits and depressions typically appear as local extrema of the mean curvature. By performing extremum search and region growing on the curvature field, the system can automatically identify these local depressions or protrusions and quantify their depth, area, and volume.

[0254] All these extracted deformation metrics will be combined with the segmented confidence assessment generated in Example 11 to form a comprehensive and traceable inspection report. As an optional approach, the extraction of deformation metrics may also include analysis of the pipe axis curvature to assess uneven settlement or displacement of the pipe.

[0255] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A three-dimensional deformable quantitative inspection method for underground pipelines, characterized in that, include: Process the raw sequences of multi-source sensors and external mapping benchmark information to generate a baseline sensor data package containing synchronous sensor sequences, wellhead event frame sets, coordinate system and external parameter metadata; The baseline sensor data packet is subjected to a first-level online correction to correct for media disturbances, and the geometrically optimized media-corrected cross-sectional data is output. By combining the baseline sensor data packet and the medium-corrected cross-sectional data, a second-level scale robustness process is performed to output the global reference pose trajectory. Based on the global reference pose trajectory and the cross-sectional data after medium correction, a three-dimensional reconstruction is performed, and a three-dimensional model and deformation quantification index are output.

2. The method according to claim 1, characterized in that, The baseline sensor data packet undergoes a first-stage online correction to correct for media disturbances, outputting geometrically optimized, media-corrected cross-sectional data, including: Synchronous sensing sequences are extracted from baseline sensing data packets, and micro-excitation sampling is performed on them using a combination of micro-amplitude angle swing and dual-power laser emission to obtain a medium state proxy quantity characterizing the medium state. A forward model of the range image with endogenous refraction and scattering correction terms is constructed. Then, using the medium state surrogate, synchronous sensing sequence, coordinate system and extrinsic data, the correction terms and cross-sectional geometric parameters are jointly solved to generate geometrically optimized medium-corrected cross-sectional data.

3. The method according to claim 2, characterized in that, The joint solution of the correction term and the cross-sectional geometric parameters includes: By utilizing synchronous sensing sequences acquired at different micro-oscillation angles and different transmission powers, a differential observation is constructed to decouple the refraction and scattering effects. Minimize the joint optimization objective function to simultaneously solve for the correction term and the cross-sectional geometry parameters. The joint optimization objective function includes: a fitting residual term that matches the difference observations, and a roundness prior penalty term based on the current cross-sectional geometry parameters.

4. The method according to claim 2, characterized in that, The method also includes: Evaluate the uncertainty of the correction term and the correction parameter of the cross-sectional geometric parameters obtained by joint solution; Based on the assessed uncertainty of the correction parameters, the micro-excitation configuration file for subsequent sampling cycles is adaptively updated. The micro-excitation configuration file controls the angular swing amplitude, power ratio, or sampling window length of the micro-excitation sampling.

5. The method according to claim 2, characterized in that, The method also includes: Multiple medium-corrected cross-sectional data generated along the pipeline mileage are used to construct an equivalent medium thickness distribution that continuously varies along the pipeline mileage and circumferential direction by introducing a joint optimization process that incorporates spatial continuity and geometric consistency priors.

6. The method according to claim 5, characterized in that, The method also includes performing at least one of the following feedback adjustment operations based on the equivalent medium thickness distribution: Feedback modulation is used for micro-excitation configuration profiles in subsequent sampling periods to enhance sampling in areas where the medium has a significant impact; When performing scale robustness processing, the weights of geometric observations are dynamically adjusted to reduce the contribution of data originating from the medium-affected region to the scale solution.

7. The method according to claim 1, characterized in that, The second-level scaling robustness is achieved by constructing and solving a unified optimization framework that integrates at least two of the following types of constraints: The dynamic consistency constraint is formed by the driving torque derived from the motor current sequence and the attachment prior. The geometric conservation constraint is extracted based on the local cross-sectional geometric stability, and the weight of the constraint is dynamically modulated by the equivalent medium thickness distribution.

8. The method according to claim 7, characterized in that, The construction of dynamic consistency constraints includes: The residual difference between the quantized wheel speed integral mileage and the inertial integral mileage; The motor current sequence is mapped to driving torque using a physical model of the motor. Establish a target term for penalizing the differential residual, where the weight of the penalty is dynamically adjusted by the mapped driving torque.

9. The method according to claim 7, characterized in that, The construction of geometric conservation constraints and their weight modulation include: Geometrically stable sections are identified from the medium-corrected cross-sectional data along the mileage to form candidates for geometric conservation constraints; The equivalent medium thickness distribution is queried. For each candidate geometric conservation constraint, a dynamic weight is assigned based on the degree of medium influence in its segment. This weight is negatively correlated with the degree of medium influence.

10. The method according to claim 1, characterized in that, The method also includes an on-demand trajectory enhancement step, including: Output a report evaluating the quality and uncertainty of the global reference pose trajectory, and automatically identify weakly observable sections based on the report; When a weakly observable section is identified, the device is triggered to perform a short-distance reverse traverse covering that section, generating bidirectional consistent observation constraints. The newly generated bidirectional consistent observation constraints are added to the unified optimization framework of scale robustness processing for secondary optimization, and the updated global reference pose trajectory is output.

Images