Method for underwater detection of deformation of structural joints of pressure tunnel

Abstract

The application discloses an underwater detection method for deformation of structural joints of a pressure tunnel structure, and comprises the following steps: controlling an underwater robot to travel along the tunnel, and acquiring real-time collected cross-section sonar data and positioning data; based on the cross-section sonar data and the positioning data, spatial registration is performed in combination with a pre-stored tunnel reference model to establish a spatial corresponding relationship between a measured cross section and a reference cross section; based on the spatial corresponding relationship, parameterized profile features of a structural joint region are extracted; and based on the parameterized profile features, deformation parameters of the structural joint are calculated, wherein the deformation parameters include a fault displacement and an opening amount. The application effectively solves the problems of large positioning cumulative error, ambiguous rotational freedom and difficulty in accurately quantifying slight deformation in underwater long-distance tunnel detection, and realizes high-precision and high-efficiency detection of structural joint deformation.

Description

Technical Field

[0001] This invention belongs to the field of intelligent sensing and monitoring of underground space engineering, and in particular, it is an underwater detection method for the deformation of structural joints in pressurized tunnels. Background Technology

[0002] Structural joints (such as expansion joints and construction joints) are weak points in the lining structure of hydraulic tunnels. Under long-term pressurized operation and geological processes, they are prone to deformations such as misalignment and opening. These minute deformations are often precursors to water-stopping failure, water jetting into the tunnel, or even structural instability. Therefore, achieving non-contact, precise, and quantitative detection of structural joint deformation without emptying the tunnel is of significant engineering and technical value for timely detection of hidden dangers and ensuring long-term structural safety.

[0003] Currently, the inspection of pressurized tunnels mainly employs remotely operated vehicles (ROVs) equipped with multibeam or cross-sectional sonar for inspection. Existing mainstream technologies typically control the ROV to navigate along the tunnel axis at a preset fixed speed, continuously scanning the inner walls with sonar to acquire point cloud data, while relying on the ROV's built-in inertial navigation system (IMU) or Doppler velocimeter (DVL) to record its trajectory. In the data processing stage, the commonly used Iterative Closest Point (ICP) algorithm is typically used to geometrically match the acquired measured point cloud with the design model. The Euclidean distance between the point clouds is calculated to generate a full-section deviation distribution map, thereby determining whether defects exist inside the tunnel.

[0004] However, the existing solutions mentioned above have technical limitations in the detection of tunnels with long distances and single geometric features, mainly manifested in three aspects: unstable positioning and registration, mismatch of sampling features, and insufficient quantization accuracy. Specifically, these include: long-distance positioning drift and ambiguity of rotational degrees of freedom. Since tunnel cross-sections are usually circular or horseshoe-shaped with high axisymmetry, and the underwater environment lacks visual texture, the general ICP registration algorithm is difficult to constrain the rotational degrees of freedom around the axis, causing the registered model to easily roll or slip along the axis. After the cumulative error of the inertial navigation system is superimposed, it is impossible to accurately restore the true spatial position of the structural joint. There is also a contradiction between the uniform sampling strategy and local small features. The width of the structural joint is usually only on the order of centimeters. If the existing uniform inspection mode is too fast, it will lead to sparse scanning cross-sections in the joint area and loss of deformation features; if the entire line is scanned at a low speed, it will generate massive amounts of redundant data and reduce work efficiency. In addition, there is the problem of deformation decoupling quantization under acoustic noise interference. Underwater sonar data itself contains a lot of reverberant noise. Simple comparison based on the original discrete point cloud not only has poor noise resistance, but also makes it difficult to accurately decouple and distinguish the two completely different physical deformation mechanisms of radial misalignment and axial opening from disordered deviation values. Summary of the Invention

[0005] The purpose of this invention is to provide an underwater detection method for the deformation of joints in pressurized tunnel structures, so as to solve the above-mentioned problems existing in the prior art.

[0006] Technical solutions, including underwater detection methods for deformation of pressurized tunnel structural joints, include:

[0007] Control the underwater robot to move along the tunnel and acquire real-time cross-sectional sonar data and positioning data;

[0008] Based on cross-sectional sonar data and positioning data, spatial registration is performed using a pre-stored tunnel reference model to establish a spatial correspondence between the measured cross-section and the reference cross-section.

[0009] Based on spatial correspondence, parametric contour features of the structural seam region are extracted;

[0010] The deformation parameters of the structural joint are calculated based on parametric contour features. The deformation parameters include misalignment and opening.

[0011] Beneficial effects: This invention effectively solves the problems of large cumulative positioning error, fuzzy rotational degrees of freedom, and difficulty in accurately quantifying minute deformations in underwater long-distance tunnel inspection, and achieves high-precision and high-efficiency detection of structural joint deformation. Attached Figure Description

[0012] Figure 1 A flowchart illustrating the steps of an underwater detection method for joint deformation in a pressurized tunnel structure, provided in an embodiment of this application.

[0013] Figure 2 A flowchart illustrating the steps for dynamically adjusting the travel speed and sonar sampling rate of an underwater robot, as provided in this embodiment of the application.

[0014] Figure 3 A flowchart illustrating the steps for establishing the spatial correspondence between the measured cross section and the reference cross section, as provided in the embodiments of this application.

[0015] Figure 4 A flowchart illustrating the steps for extracting parametric contour features of a structural seam region, as provided in an embodiment of this application. Detailed Implementation

[0016] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0017] It should be noted that the terms include and have, and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0018] like Figure 1 As shown, an underwater method for detecting the deformation of joints in pressurized tunnel structures includes the following steps:

[0019] Control the underwater robot to move along the tunnel and acquire real-time cross-sectional sonar data and positioning data.

[0020] In this embodiment, the underwater robot (ROV) can be a medium or small-sized operational remotely operated vehicle with current resistance and stable attitude control capabilities, and it must be equipped with a high-precision cross-sectional sonar system. For example, a single-beam mechanical scanning sonar or a multi-beam imaging sonar can be used, mounted on the front or top of the ROV, with the scanning plane perpendicular to the tunnel axis, to acquire distance information reflecting the contour of the tunnel's inner wall. The cross-sectional sonar data is typically output in polar coordinates (distance ρ, angle θ), which can then be converted into point cloud data in a Cartesian coordinate system.

[0021] To achieve precise positioning, ROVs also need to be equipped with a combined navigation system. Positioning data is specifically calculated by fusing data from a Doppler velocity level (DVL), an inertial measurement unit (IMU), and a depth gauge. The DVL provides bottom velocity, the IMU provides three-axis acceleration and angular velocity, and the depth gauge provides high-precision depth information. The positioning data is specifically represented as a position vector S including a timestamp. pos (t)=[x(t),y(t),z(t),φ(t),θ(t),ψ(t)] T Where x, y, and z represent the east, north, and sky positions, respectively; φ, θ, and ψ represent the roll, pitch, and yaw angles, respectively; and t represents the time variable. T This is a transpose.

[0022] In some alternative implementations, considering the complexity of the underwater environment, the system can automatically switch positioning modes if a single sensor fails. For example, when the DVL loses its bottom lock, visual odometry-assisted positioning can be temporarily performed using sonar image feature matching, or coarse correction can be made using cable length counter data to ensure the continuity of positioning data. This multi-source fusion positioning method provides reliable initial position information for subsequent spatial registration.

[0023] Based on cross-sectional sonar data and positioning data, spatial registration is performed using a pre-stored tunnel reference model to establish a spatial correspondence between the measured cross-section and the reference cross-section.

[0024] Specifically, the measured data collected by the ROV, which contains noise and attitude deviations, is aligned with the ideal design model to the same coordinate system. The tunnel reference model refers to a pre-constructed three-dimensional digital model that reflects the geometric shape of the tunnel design, such as a standard cylindrical or horseshoe-shaped point cloud model derived from BIM (Building Information Modeling). Spatial registration is not merely a simple coordinate transformation, but a process of solving for the optimal rigid body transformation matrix, namely the rotation matrix R and the translation vector t', to maximize the spatial overlap between the measured point cloud and the reference model.

[0025] Because underwater sonar data is typically collected in the ROV's own volume coordinate system and is affected by water flow, ROVs inevitably experience yaw or roll during their movement. Without spatial registration, directly analyzing point cloud data can easily misinterpret the ROV's tilt as deformation of the tunnel cross-section. By establishing a spatial correspondence, all measured data can be uniformly mapped to the tunnel's global coordinate system, eliminating observation errors caused by ROV movement.

[0026] In some alternative implementations, to improve the efficiency and success rate of registration, coarse registration can be performed first using approximate locations in the positioning data to roughly move the measured data to the vicinity of the reference model; then, fine registration algorithms such as the Iterative Closest Point (ICP) algorithm and its variants can be used for fine-tuning. Furthermore, for long-distance tunnels with simple features, conventional registration is prone to errors such as axial slippage or rotation around the axis. Therefore, a segmented registration strategy based on geometric constraints is preferred, utilizing the unique floor and structural joint features of the tunnel to lock the degrees of freedom.

[0027] Based on spatial correspondence, parametric contour features of the structural seam region are extracted.

[0028] In this embodiment, the structural seam area is accurately extracted from massive point cloud data and converted into mathematically computable parametric form. After spatial registration, the measured data is aligned with the benchmark model. At this point, based on the pre-established structural seam design mileage information, point cloud data of the structural seam and a certain range before and after it can be extracted from the entire point cloud. This area is the structural seam region. Extracting parametric contour features means no longer directly using discrete, noisy raw point clouds for subsequent calculations, but constructing a continuous curve or function describing the cross-sectional shape through mathematical fitting methods. For example, the least squares method can be used to fit a circle or an ellipse to the cross-sectional point cloud. Preferably, in order to accurately describe the local deformation at the structural seam, such as misalignment or spalling, spline curves (such as B-splines) can be used for fitting. The parametric contour features are specifically represented by the control point sequence P of the fitted curve. ctrl =[P0, P1, ..., P mThese control points not only preserve the macroscopic shape information of the cross-section, but also suppress sonar measurement noise through smoothing, providing a stable mathematical basis for the quantification of minute deformations.

[0029] In some alternative implementations, the feature extraction strategy can be flexibly adjusted for tunnel cross-sections of different shapes. For example, for a horseshoe-shaped cross-section, segmental fitting can be performed on the floor, sidewalls, and crown separately; for a circular cross-section, overall closed curve fitting can be performed. Regardless of the fitting method used, the core objective is to transform unstructured point clouds into structured geometric parameters.

[0030] The deformation parameters of the structural joint are calculated based on parametric contour features. The deformation parameters include misalignment and opening.

[0031] Specifically, by comparing the parametric profile features of the cross-sections before and after the structural joint, such as comparing the coordinates of corresponding control points on both sides, the specific deformation can be analyzed. Deformation parameters typically include the misalignment amount d. offset The difference in height and the opening d between the two sides of the seam in the radial direction. open That is, the change in axial distance between the two sides of the seam and the possible relative rotation θ rot .

[0032] For example, if the control points on the side behind the seam are found to be 5 mm higher in the radial direction than those on the side in front of the seam, it can be determined that there is a 5 mm misalignment deformation at that location. The calculation method based on parametric features avoids misjudgments caused by uneven point cloud density or noise interference during direct point-to-point comparisons. The calculated deformation parameter D=[d offset d open θ rot ] T This will serve as the final test result, used to assess the health status of the tunnel structure.

[0033] In some optional implementations, the system can also generate visual analysis charts based on the calculated deformation parameters, such as the distribution curve of misalignment along the tunnel mileage and the contour comparison diagram of typical defect sections, to intuitively assist engineers in making decisions. Simultaneously, if the deformation parameters exceed a preset safety threshold, the system can automatically trigger an early warning mechanism.

[0034] This embodiment solves the problems of large cumulative positioning error, high rotation uncertainty, and difficulty in quantifying deformation in long-distance tunnel detection in the prior art by autonomously collecting data in the data acquisition stage and introducing spatial registration and parametric feature extraction in the data processing stage, thus realizing accurate and safe detection of structural joint deformation.

[0035] In one possible implementation, before acquiring the real-time cross-sectional sonar data and positioning data, a structural joint mileage distribution database is pre-built, specifically as follows:

[0036] Obtain the as-built drawings or design data of the tunnel and extract the design mileage locations of all structural joints.

[0037] In this embodiment, the as-built drawings or design documents typically include longitudinal sections of the tunnel, pipe section layout diagrams, and detailed joint design drawings. Workers can use digital means to access this data and extract the mileage marker for each structural joint in the design coordinate system. The design mileage location is denoted as L. design It is usually expressed in meters, representing the axial distance of the structural joint from the tunnel entrance or a specific reference point.

[0038] In the specific implementation process, if the drawings are electronic CAD files, scripts can be written to automatically search for specific layers or block attributes and extract the coordinate information representing joints. If the drawings are scanned copies of paper documents, optical character recognition (OCR) technology combined with manual verification can be used to identify mileage markings on the drawings. For example, for a 1000-meter-long water conveyance tunnel, the design drawings may show a construction joint every 10 meters and an expansion joint every 50 meters. By extracting this information, a preliminary list of theoretical locations can be constructed.

[0039] In some alternative implementations, considering the possible changes during tunnel construction, if construction logs or supervision records can be obtained, the extracted design mileage location can be corrected based on these auxiliary materials to ensure the authenticity of the data.

[0040] Alternatively, obtain historical inspection data of the tunnel and identify the measured mileage location of structural joints in the historical data.

[0041] Specifically, it provides a more accurate data source, particularly for older tunnels that have been in operation for many years. Historical inspection data can include past manual inspection records, historical sonar scan point clouds, or underwater camera videos. By analyzing this data, the physical location of the structural joint in the actual environment can be directly observed, thereby obtaining the measured mileage location, denoted as L. measured .

[0042] For example, in historical sonar point cloud data, structural joints typically appear as circumferential grooves or protrusions on the pipe wall. By replaying historical point clouds and identifying features, the mileage of these geometric features can be marked. Because this data comes directly from the constructed entity, it includes all construction errors and positional variations caused by geological settlement, and its accuracy is generally superior to design drawings.

[0043] In some alternative implementations, if the historical data includes visual images, image recognition algorithms can be used to detect the texture features of cracks or joints, and the mileage of the structural joint can be calculated by combining this with the location records at the time. For tunnels lacking digitized historical data, station information from manual maintenance logbooks can also be used.

[0044] The designed or measured mileage locations are compiled to generate a structural joint mileage distribution database, which serves as the source for the structural joint location list.

[0045] Specifically, the multi-source data is cleaned, fused, and structured. To ensure database reliability, the system needs to formulate a data fusion strategy. Generally, if the same structural joint has both designed mileage and historical measured mileage, the measured mileage location L is preferred. measured This is the final record because it better reflects the current state of the tunnel; if only design data is available, then the design mileage location L should be used. design The structural joint mileage distribution database can be stored using relational database tables or lookup tables. Each record in the table corresponds to one structural joint, and the key fields include at least the structural joint's unique identifier ID. seam Mileage L And optional joint type, such as construction joint or expansion joint. For example, a record in the database can be represented as: {ID: 001, Mileage: 105.5, Type: Expansion}, where Expansion represents an expansion joint.

[0046] In some optional implementations, the generated database may also contain prior information on the geometric properties of the structural joints, such as the design joint width or pipe section diameter, for subsequent algorithmic use. The final generated database will be downloaded to the underwater robot's control computer as a list of structural joint locations, used to calculate the real-time joint spacing d during real-time detection. seam (t) guides the robot's speed control and sonar sampling strategies. The offline construction and online invocation model reduces the complexity of real-time computation while ensuring the relevance of the detection strategy.

[0047] This embodiment serves as a preliminary preparation step, providing necessary prior data support for the adaptive acquisition based on real-time seam spacing and the segmented registration based on geometric constraints in subsequent embodiments. By integrating design data and historical measured data, it can effectively solve the problem of inaccurate positioning benchmarks caused by construction deviations or missing drawings.

[0048] like Figure 2 As shown, in an exemplary embodiment, after acquiring the real-time collected positioning data, the method further includes dynamically adjusting the underwater robot's travel speed and sonar sampling rate, specifically:

[0049] Based on positioning data, the real-time gap between the underwater robot's current position and the nearest structural gap in the pre-stored list of structural gap locations is calculated.

[0050] In this embodiment, the system needs to acquire high-precision real-time mileage L(t). Although the positioning data contains location information, the mileage obtained by direct integration may be inaccurate due to the cumulative error of the underwater Doppler velocity meter (DVL). Therefore, it is preferable to use a multi-source information fusion method to calculate the current mileage. Specifically, the real-time mileage is estimated using the following formula:

[0051] L(t) = α' * L DVL (t)+(1-α')*L feature (t);

[0052] Where L(t) is the current mileage after fusion, L DVL (t) represents the mileage obtained by integrating the velocity measurements from a Doppler velocimeter, L feature (t) represents the corrected mileage obtained based on sonar image features or visual odometry matching, and α' is the fusion weighting coefficient, typically ranging from 0 to 1, depending on the real-time confidence level of each sensor. Specifically, α' can be dynamically calculated based on the ratio of the variance of DVL velocity measurement to the variance of feature-matched mileage, or a fixed value can be preset by comparing data from calibrated flight segments before actual operation. Those skilled in the art can determine appropriate values ​​based on the actual sensor performance and tunnel environmental conditions.

[0053] After obtaining the current mileage L(t), the system calls the pre-built structural joint mileage distribution database and traverses all structural joint mileage locations L in the database. seam_k It calculates the absolute distance between the current position and the nearest structural seam in the list. Real-time seam distance d seam The formula for calculating (t) is as follows:

[0054] d seam (t)=min(abs(L seam_k -L(t)));

[0055] Here, min represents the minimum value function, and abs represents the absolute value function. By calculating the real-time seam distance, the system can sense in real time how far away the next structural seam that needs to be focused on for inspection is.

[0056] The underwater robot's travel speed and sonar sampling rate are dynamically adjusted based on the real-time gap. When the real-time gap is less than the preset encryption threshold, the travel speed is reduced and the sonar sampling rate is increased.

[0057] Specifically, the encryption threshold W analysisThe range of the core concern area of ​​the structural joint is defined, for example, it can be set to 0.5 meters. When the real-time joint spacing d... seam When (t) exceeds this threshold, the system determines that the underwater robot is in a normal pipe segment. In this case, it maintains a high traveling speed to improve operational efficiency and maintains a standard sonar sampling rate to save storage space. When the real-time slot distance d... seam When (t) is less than or equal to the encryption threshold, the system determines that it has entered the key detection zone of the structural seam. At this time, the control system automatically issues a deceleration command to reduce the thruster power of the underwater robot, and at the same time sends a command to the cross-section sonar to increase the trigger frequency. This increases the number of cross-section scans per unit length. For example, let the travel speed be v and the sampling rate be f. s Then the sampling interval along the tunnel axis ΔL*=v / f s By simultaneously decreasing v and increasing f s This can reduce the sampling interval ΔL*, thereby obtaining a sufficient number of cross-sectional data within a limited seam width range, providing ample samples for subsequent deformation calculation.

[0058] In some alternative implementations, dynamic adjustment is not a step-like on / off control, because sudden changes in speed would cause drastic fluctuations in the underwater robot's pitch angle, resulting in wobbling of the sonar scanning plane. Therefore, a continuously varying function is preferred to achieve a smooth transition.

[0059] In a preferred implementation, the underwater robot's travel speed and sonar sampling rate are dynamically adjusted based on the real-time slot gap. Specifically, the target travel speed and target sampling rate of the underwater robot are calculated using the following formula:

[0060] v(t) = v max -(v max -v min )*exp(-(d seam (t)) 2 / (2*σ v 2 ));

[0061] f s (t)=f min +(f max -f min )*exp(-(d seam (t)) 2 / (2*σ f 2 ));

[0062] Where v(t) is the target speed, f s (t) represents the target sampling rate, d seam (t) represents the real-time seam distance, v max and v minThese are the preset maximum and minimum speeds, f. max and f min These are the preset maximum and minimum sampling rates, σ v and σ f The preset adjustment range parameter is exp(), which is an exponential function.

[0063] In this embodiment, a specific mathematical model for adaptive adjustment is presented, employing a variant of the Gaussian function to achieve smooth parameter decay and recovery. For example, v max It can be set to 0.5 meters per second, v min Set to 0.1 meters per second. max It can be set to 10 Hz, f min Set to 2 Hz. σ v and σ f The sensitivity range used to control changes in speed and sampling rate is recommended to be 5 to 10 times the width of the structural joint, for example, set to 0.2 meters.

[0064] To illustrate this more intuitively, let's take speed control as an example and explain the numerical values: When the underwater robot is directly facing the center of the structural seam, the real-time seam distance d... seam When (t) is 0, the exponent term exp(0) equals 1, and the calculated result of the target velocity v(t) is v. min That is, at a speed of 0.1 meters per second, reaching the slowest cruising state, ensuring a precise scan of the gap. When the underwater robot moves away from the structural gap, for example, d... seam (t) is much greater than σ v When the exponential term approaches 0, the target velocity v(t) approaches v0. max That is, at 0.5 meters per second, rapid cruise is resumed. The control curve based on the Gaussian function has good smoothness, avoiding impact on the mechanical structure.

[0065] In a further embodiment, dynamically adjusting the underwater robot's travel speed also includes applying acceleration limiting processing to the target travel speed, specifically:

[0066] Calculate the rate of change of the target's speed between adjacent time points;

[0067] When the absolute value of the rate of change exceeds the preset acceleration threshold, the target travel speed is limited by filtering to keep the rate of change of the underwater robot's speed within the acceleration threshold range, so as to maintain the attitude stability during sonar acquisition.

[0068] In this embodiment, although the Gaussian function described above is theoretically smooth, in practical engineering, jumps in positioning data or delays in the control cycle may cause sudden changes in the calculated target velocity. To prevent pitch oscillations caused by rapid acceleration or deceleration of the underwater robot, dynamic constraints are required on the velocity command.

[0069] Specifically, the system calculates the rate of change 'a' of the target's speed between adjacent time points. cmd The calculation formula is as follows:

[0070] a cmd =(v target (t)-v cmd (t prev )) / Δt;

[0071] Among them, v target (t) represents the theoretical target velocity at the current moment, v cmd (t prev ) represents the actual command speed issued in the previous control cycle, and Δt represents the time interval of the control cycle.

[0072] The rate of change is compared with a preset acceleration threshold a. max Compare the acceleration threshold a. max Based on the dynamic characteristics of the underwater robot, for example, a threshold of 0.05 meters per second squared. If the absolute value of the calculated rate of change exceeds this threshold, a limiting filter is initiated. The specific limiting filter logic is as follows:

[0073] v cmd (t)=v cmd (t prev )+sign(a cmd )*min(abs(v target (t)-v cmd (t prev )), a max *Δt);

[0074] Where sign is the sign function, min is the minimum value function, and v cmd (t) represents the actual speed of instructions issued during the current control cycle.

[0075] The above processing ensures that the command speed v sent to the underlying thruster is maintained. cmd The rate of change of (t) is always confined to [-a] max +a max Within the specified range, attitude disturbances during ROV acceleration and deceleration were effectively suppressed, ensuring that the cross-sectional sonar scanning plane remained perpendicular to the tunnel axis, thus improving the geometric consistency of the acquired data.

[0076] This embodiment calculates the relative position of the underwater robot to the structural seam in real time, and automatically reduces its travel speed and increases the sonar sampling frequency when approaching the structural seam. This improves the point cloud density and data quality of the key area without significantly reducing the overall inspection efficiency.

[0077] like Figure 3 As shown, according to one aspect of this application, establishing a spatial correspondence between measured cross-sections and reference cross-sections includes:

[0078] Based on the pre-stored structural joint mileage distribution, the cross-sectional sonar data is divided into predetermined independent structural segment point clouds.

[0079] In this embodiment, since the errors of the underwater robot's inertial navigation system accumulate linearly or nonlinearly over time, directly registering the entire point cloud spanning several kilometers often fails to converge. Therefore, a strategy of breaking down the problem into smaller parts can be adopted. Specifically, the system calls a pre-built structural joint mileage distribution database to obtain the design mileage L of each structural joint. seam_k Using two adjacent structural joints L seam_k and L seam_(k+1) Using this as a boundary, the continuously acquired cross-sectional sonar data is truncated on the time axis or approximate mileage axis to form independent structural segment point clouds P. k Each structural segment point cloud P k This corresponds to a specific pipe segment. Segmented processing not only isolates the cross-segment propagation of errors and eliminates the impact of long-distance cumulative errors, but also allows each pipe segment to have independent rigid body transformation parameters, thereby adapting to uneven settlement of the tunnel or misalignment between pipe segments.

[0080] For each structural segment point cloud, a registration objective function is constructed. The registration objective function includes a data fitting term to minimize the point cloud matching error and a geometric feature constraint term to limit the degrees of freedom of rigid body transformation.

[0081] Specifically, the traditional Iterative Closest Point (ICP) algorithm only includes a data fitting term, which aims to minimize the distance from a point to a surface. However, in a smooth circular tunnel, when the point cloud slides along or rotates around the axis, the distance from a point to a surface may change very little, causing the algorithm to get stuck in a local optimum. To address this, this embodiment adds geometric feature constraints, namely, the base plate plane constraint and the longitudinal position correction constraint, to the registration objective function.

[0082] Building upon this, to further improve the robustness of registration, the registration objective function also includes adaptive weight allocation based on region features, specifically:

[0083] Based on the cross-sectional geometric features, the points in the point cloud of the structural segment are classified into points in the regular pipe section area, points in the structural joint area, or points in the abnormal deformation area.

[0084] In the data fitting term of the registration objective function, different confidence weights are assigned to points of different categories; among them, the confidence weight of points in the regular pipe section area is higher than that of points in the structural joint area and points in the abnormal deformation area.

[0085] Specifically, the local curvature or normal rate of change of each point in the point cloud of the structural segment is calculated as a cross-sectional geometric feature. The system classifies the points in the structural segment point cloud based on these cross-sectional geometric features. For example, the normal vector curvature or roughness of each point is calculated, and regions with smooth curvature are marked as regular pipe section points, while regions with abrupt curvature changes are marked as structural joint points or abnormal deformation points. When constructing the data fitting term, different confidence weights w are assigned to points of different classifications. i Typically, the geometry of regular pipe sections best approximates the design model, and therefore receives a higher weight, for example, w. high =1.0; however, due to the presence of spalling, attachments, or actual deformation in structural joint areas or abnormal deformation zones, their point cloud positions may deviate from the design model. Assigning high weights to these areas would skew the overall registration results; therefore, lower weights are given, for example, w. low =0.1. Through a differentiated weighting strategy, the registration process is mainly driven by reliable pipe wall data, without being affected by local deformation.

[0086] For example, the specific implementation of the region classification and differentiated registration strategy is as follows: the system divides the point cloud of the structural segment into three types of regions according to the cross-sectional geometric features, and adopts different registration parameter settings for each type of region.

[0087] The first category is the regular pipe section region, characterized by a smooth cross-sectional profile, a geometric shape close to the design value, and small local curvature variations. For this region, strong geometric constraints are adopted in the registration objective function, namely, increasing the values ​​of the bottom plate constraint weight λ1 and the longitudinal constraint weight λ2, while using a larger iteration step size to accelerate the convergence speed.

[0088] The second category is structural seam regions, characterized by significant geometric discontinuities, including grooves, steps, or gaps. For this region, a weak geometric constraint is adopted in the registration objective function, i.e., reducing the values ​​of λ1 and λ2 to avoid misclassifying true structural seam features as registration errors and forcibly flattening them. Simultaneously, a smaller iteration step size and a more stringent convergence criterion are used to improve registration accuracy.

[0089] The third category is the abnormal deformation area, characterized by local depressions, protrusions, peeling, or attachments. For this area, a robust registration strategy is adopted, introducing the Huber loss function or truncated least squares into the data fitting term to reduce the bias of outliers on the overall registration.

[0090] The criteria for region classification can be set based on the local curvature variance or the deviation from the design model. For example, when the local curvature variance is less than a threshold σ... κ_th And the deviation is less than the threshold d th When the deviation exceeds the threshold but meets the geometric characteristics of a structural joint, it is determined to be a regular pipe section area; when the deviation exceeds the threshold but meets the geometric characteristics of a structural joint, it is determined to be a structural joint area; otherwise, it is determined to be an abnormal deformation area.

[0091] In one embodiment of this application, the geometric feature constraint includes a base plate plane constraint, the construction process of which includes:

[0092] Extract the base plate region point set from the structural segment point cloud, and use a plane fitting algorithm to calculate the measured base plate normal vector of the base plate region point set;

[0093] Obtain the corresponding design floor normal vector from the tunnel reference model;

[0094] The constraints on the bottom plate plane are constructed to minimize the angle between the measured bottom plate normal vector and the designed bottom plate normal vector, and the rotational degrees of freedom about the tunnel axis are locked.

[0095] In this embodiment, for most water conveyance tunnels, although the cross-section may be circular or horseshoe-shaped, the bottom is usually designed with a flat base or a naturally deposited horizontal surface, serving as a natural level. The system extracts point clouds within the bottom angle range of the structural segment point cloud as a candidate base region point set. To eliminate potential rocks, silt, or abnormal noise points at the bottom, the Random Sample Consensus (RANSAC) algorithm is preferably used for plane fitting. The RANSAC algorithm constructs a plane model ax + by + c * z + d = 0 by randomly sampling a subset, counts the number of local points, and selects the model with the most local points as the optimal base plane after multiple iterations, where a, b, and c are the three components of the plane normal vector, and d is the plane's position parameter.

[0096] The plane equation coefficients obtained from the fitting can be used to determine the measured base plate normal vector n. meas =[a, b, c] T Meanwhile, it can be seen from the design drawings or benchmark model that the ideal base plate normal vector n design Typically perpendicular to the geoid, for example, [0, 1, 0] T By adding a term to the objective function, n is forced to... meas With n design Parallelism, meaning that the cross product of two vectors is zero or the dot product is 1, can effectively lock the rotational degree of freedom of the point cloud around the tunnel axis, preventing the registered pipe section from becoming skewed.

[0097] In another embodiment of this application, the geometric feature constraint term further includes a longitudinal position correction constraint, the construction process of which includes:

[0098] Calculate the rate of change of profile curvature of each section in the point cloud of the structural segment, and identify the location where the rate of change of profile curvature exceeds the preset abrupt change threshold as the actual structural seam location;

[0099] Obtain the design structural joint position corresponding to the measured structural joint position, and calculate the longitudinal deviation between the two.

[0100] Constraints on the longitudinal position are constructed based on the longitudinal deviation value, the mileage component in the positioning data is corrected, and the translational degree of freedom along the tunnel axis is locked.

[0101] Specifically, because the sonar beam has a certain width, the depth value of the echo signal will abruptly change when scanning through structural seams. This abrupt change manifests as an extreme value of curvature in differential geometry. Before detecting curvature abrupt changes, axial density analysis can be used as a preliminary screening method for the location of structural seams. Axial density analysis refers to counting the number of point clouds per unit length along the tunnel axis. Due to the presence of grooves or joint gaps at structural seams, sonar echoes may show signal loss or abnormal enhancement in this area, leading to abrupt changes in point cloud density. Specifically, the tunnel axis is divided into sections at fixed intervals Δz. bin Divide the data into several statistical intervals and calculate the point cloud density ρ in each interval. density (z), the formula is as follows:

[0102] ρ density (z) = N points (z) / Δz bin ;

[0103] Where, ρ density (z) represents the axial point cloud density at mileage z, N points (z) represents the area centered at z with a width of Δz. bin The number of point clouds within the interval, Δz bin This is the width of the statistical interval.

[0104] When ρ density (z) When the density value decreases or increases significantly relative to the adjacent intervals before and after it, a structural seam may exist at that location. Using the location of the density abrupt change as a candidate structural seam location, and then combining it with the subsequent curvature change rate detection for precise localization, can improve the robustness of structural seam identification.

[0105] The system calculates the rate of curvature change κ'(z) of each cross-sectional profile along the tunnel axis z. Preferably, the calculation of the cross-sectional profile curvature uses a discrete point numerical approximation method. For the ordered point sequence formed by sorting the cross-sectional point cloud according to polar coordinate angles, three adjacent points p are taken. (i-1) p i p (i+1) The arc formed by this point p has a curvature approximately equal to that of the point p. i Local curvature κ(θ) at i Specifically, using the geometric principle that a circle is determined by three points, point p... i The formula for calculating the curvature at a given point is as follows:

[0106] κ(θ i ) = (4 * A triangle) / (L a * L b * L c );

[0107] Wherein, κ(θ) i ) represents the polar coordinate angle θ i The curvature value at point A triangle For three points p (i-1) p i p (i+1) The area of ​​the enclosed triangle, L a For point p (i-1) With p i The distance between them, L b For point p i With p (i+1) The distance between them, L c For point p (i-1) With p (i+1) The distance between them.

[0108] After obtaining the curvature sequence distributed circumferentially along the cross-section, the rate of change of curvature along the axial direction z, κ'(z), is further calculated. For a cross-section located at mileage z, the curvature difference between this cross-section and its adjacent cross-sections at the same angular position is taken to obtain the rate of change of curvature at that location. When the absolute value of κ'(z) exceeds a preset abrupt change threshold and exhibits a zero-crossing characteristic of first positive then negative or first negative then positive, the location is determined to be the measured structural joint location L. meas The preset mutation threshold can be determined by calibrating historical data at known structural joint locations based on the actual tunnel cross-sectional curvature characteristics.

[0109] As an alternative approach, structural joints can also be identified using the shape difference index ΔS(z) between adjacent sections. The calculation formula is as follows:

[0110] ΔS(z)=(1 / N)*Σ||p i (z)-p i (z-Δz)||;

[0111] Where N is the number of cross-sectional points, p i (z) represents a point on the current cross-section, p i (z-Δz) represents the corresponding point on the previous cross-section, and Σ represents the summation of i from 1 to N. When the sonar sweeps across the structural joint, the cross-sectional shape changes drastically, and ΔS(z) will show a significant local peak. The mileage corresponding to this peak is the measured location L of the structural joint. meas .

[0112] Regardless of whether the curvature method or the shape difference method is used, once the measured structural joint location L is determined... meas Compare it with the design mileage L in the database designBy comparing the values, the longitudinal deviation ΔL = L can be calculated. design -L meas The longitudinal deviation value represents the cumulative error of the inertial navigation system at that location. By introducing deviation constraints into the objective function, it is equivalent to using structural seams to adjust the sliding point cloud to the correct design position.

[0113] By iteratively optimizing the registration objective function, the spatial transformation matrix that minimizes the registration objective function is obtained, which serves as the spatial correspondence.

[0114] In this embodiment, the registration objective function E is a multi-objective weighted function. In a preferred implementation, the registration objective function is solved by iterative optimization, specifically by finding the minimum value of the following registration objective function to determine the rotation matrix R and the translation vector t':

[0115] E=Σ||R·p i +t'-q i || 2 +λ1·||R·n meas -n design || 2 +λ2·||t z -ΔL|| 2 ;

[0116] Where E is the registration objective function value, p i For points in the structural segment point cloud, q i For the tunnel reference model and p i The corresponding nearest point, n meas Let n be the measured base plate normal vector. design To design the base plate normal vector, t z Let λ' be the component of the translation vector t' along the tunnel axis, ΔL be the longitudinal deviation value, λ1 and λ2 be preset constraint weight coefficients, and ||| represent the Euclidean norm. In a preferred embodiment, λ1 is set to 10.0 and λ2 is set to 5.0. The specific values ​​of λ1 and λ2 can be adaptively adjusted according to the tunnel cross-sectional geometry, the reliability of the floor plan fitting, and the accuracy of structural joint identification, for example, by optimizing parameters on a verification section with known registration results.

[0117] In another possible embodiment, the registration objective function uses a weighted point-to-surface distance metric, specifically expressed as follows:

[0118] E=Σ(w i *((R*p i +t'-q i )·n i ) 2 )+λ1*norm(R*n meas -ndesign ) 2 +λ2*(t z -ΔL) 2 ;

[0119] Where E is the total energy value of the registration objective function, Σ represents the summation over all points in the structural segment, and w i Let p be the adaptive weight of the i-th point, R be the 3x3 rotation matrix to be determined, and p be the adaptive weight of the i-th point. i Let be the coordinate vector of the i-th point in the point cloud of the structural segment, t' be the three-dimensional translation vector to be determined, and q be the coordinate vector of the i-th point in the point cloud. i For the reference model and the transformed point p i The coordinate vector of the corresponding nearest neighbor, n i For the nearest neighbor point q i The unit normal vector at n, the symbol · represents the vector dot product operation, n meas Let n be the measured base plate normal vector. design To design the base plate normal vector, norm represents the Euclidean norm of the vector, t z Let t' be the component of the translation vector along the tunnel axis, ΔL be the longitudinal deviation correction value, and λ1 and λ2 be the preset constraint weight coefficients.

[0120] Compared to point-to-point distance metrics, point-to-surface distance metrics only penalize deviations along the normal direction, allowing points to slide in the tangential direction, thus exhibiting better convergence for smooth surfaces. The registration algorithm employs an iterative solution strategy, with each iteration including nearest neighbor search and transformation parameter update, until the change in the objective function is less than a preset threshold or the maximum number of iterations is reached.

[0121] Optionally, the registration objective function is a nonlinear least squares problem, which can be solved iteratively using the Levenberg-Marquardt algorithm or the Gauss-Newton method. The final output optimal transformation parameters T*=[R*, t'*] represent the spatial correspondence of the structural segment. By iteratively optimizing the registration objective function, each point cloud segment is precisely placed in its true physical location with a correct orientation, laying a solid geometric foundation for subsequent extraction of minor deformations.

[0122] This embodiment proposes to transform the global registration problem into a piecewise local optimization problem, and introduces base plate normal constraints and longitudinal position correction, upgrading the traditional pure data-driven registration to a dual-driven registration of geometry and data, thereby improving the physical authenticity of the registration results.

[0123] like Figure 4 As shown, in one embodiment of this application, extracting the parametric contour features of the structural seam region includes:

[0124] By utilizing spatial correspondence, the cross-sectional sonar data are unified to the tunnel reference coordinate system, and the point cloud of the structural joint area is sliced ​​along the tunnel axis to obtain a predetermined set of cross-sectional points.

[0125] Alternatively, the spatial correspondence is used to unify the cross-sectional sonar data to the tunnel reference coordinate system, and the structural joint area is determined according to the pre-stored structural joint mileage location. The point cloud of the structural joint area is sliced ​​along the tunnel axis to obtain a predetermined set of cross-sectional points.

[0126] Specifically, after fine registration, all sonar point cloud data has been transformed into a unified tunnel reference coordinate system, whose Z-axis typically coincides with the tunnel's design axis. The system then uses the mileage L of the structural joint as the reference coordinate system. seam Define the analysis window, for example, select [L] seam -0.5, L seam Within a range of +0.5 meters, all point cloud data within this range are extracted. Further, to analyze the axial variation of the cross-sectional shape, the 3D point cloud needs to be converted into 2D slices. The system sets a fixed slice spacing Δz, for example, 0.05 meters, and extracts a thin slice of point cloud with a thickness of δz every Δz along the Z-axis. The points within each slice are projected onto the corresponding XY plane, forming several independent 2D cross-sectional point sets P. section_k Slicing simplifies complex spatial surface analysis into a series of planar curve analyses, reducing computational dimensionality and making it easier to capture abrupt changes along the axis.

[0127] For each cross-sectional point set, B-spline curve fitting is performed, and a pre-defined uniform node vector configuration is used during the fitting process.

[0128] In this embodiment, traditional circle fitting can only obtain the radius and center, failing to describe local concavities or non-circular deformations; while directly connecting discrete points retains all sonar measurement noise. Therefore, it is preferable to use a cubic uniform B-spline curve C(u) to approximate the cross-sectional point set. The B-spline curve consists of a set of control points P j and basis functions N j，k (u) is defined as follows: C(u) = Σ(P j *N j，k (u)); where Σ represents the summation over control point index j, and k is the curve order, which is 4 for cubic B-spline curves.

[0129] To ensure the comparability of fitted curves at different locations or at different times, a uniform node vector configuration must be enforced. The node vector U = [u0, u1, ..., u...]. m+k+1The distribution of the basis functions is determined. For example, the node vectors are fixed to a uniform distribution, and the number of nodes is predetermined based on the cross-sectional complexity. For instance, for a standard circular cross-section, 8 to 12 control points are preferably set; for a horseshoe-shaped or archway-shaped cross-section, 16 to 24 control points are preferably set. This uniform configuration ensures that the control points P... j They have a clear correspondence in geometric position; for example, P0 always corresponds to the top of the section, and P4 always corresponds to the right wall.

[0130] In one possible implementation, B-spline curve fitting is performed on each cross-sectional point set, specifically by solving an optimization objective function that includes data approximation terms and smoothness regularization terms:

[0131] minΣ||p l -C(u l )|| 2 +λ smooth ·Σ||P (j+1) -2P j +P (j-1) || 2 ;

[0132] Where min represents the minimization operation, p l For the measured points in the cross-section point set, C(u) l ) represents the corresponding point on the B-spline curve, used to ensure that the curve is as close as possible to the measurement data; P j Let λ be the control point. smooth The smoothness weighting coefficient is a preset value. The smoothness regularization term is used to suppress local oscillations in the curve caused by sonar noise.

[0133] Specifically, underwater sonar data often exhibits rough point clouds due to multipath effects and reverberation noise. To prevent overfitting (non-physical oscillations) in the fitted curve, a second-order difference regularization constraint is introduced when solving for the control points. The system solves for the minimum of the optimization objective function, where the first term is the data approximation term. The second term is the smoothness regularization term, which uses the second derivative (curvature) of the second-order difference approximation curve to penalize severe jitter in the control point sequence. The smoothness weight coefficient λ... smooth The larger the value, the smoother the fitted curve and the stronger its noise resistance. By solving this linear least squares problem, a set of control point sequences P that closely follows the measurement data and has good smoothness can be obtained. j .

[0134] In one possible implementation, the basis functions N in the B-spline curve j，k (u) is a set of piecewise polynomial functions defined on the parameter interval. It has local support, meaning that each basis function is non-zero only in a few adjacent node intervals, such that the control point P jIt only affects the local shape of the curve. For example, the B-spline basis functions are defined using the Cox-deBoor recursive formula, employing uniformly periodic node vectors for closed cross-sectional profiles.

[0135] This embodiment uses a cubic B-spline curve, where k is 4, and the corresponding basis function is a piecewise cubic polynomial, exhibiting second-order continuous differentiability and smoothness. The node vectors are uniformly distributed, meaning the spacing between adjacent nodes is equal, ensuring that the influence of each control point on the curve shape is fairly balanced, facilitating subsequent comparisons of data from multiple periods.

[0136] It should be noted that the design of the node vector must satisfy the following constraints: for a k-order B-spline curve with m+1 control points, the length of the node vector must be m+k+1. The node values ​​in the node vector determine the distribution of the support intervals of each basis function, thus affecting the parameterization method of the curve. For closed cross-sectional profile curves, a uniform periodic node vector can be used. The specific construction method is as follows: divide the parameter interval [0, 1] into m+1 sub-intervals, and take the node values ​​as i / (m+1) sequentially, where i traverses from 0 to m+k. Repeat k-1 nodes at the beginning and end of the curve to achieve closed connection. The selection principle for the number of control points m+1 is: for geometrically simple circular cross-sections, m of 8 to 12 is sufficient to fully express the cross-sectional shape; for geometrically complex horseshoe or archway cross-sections, m of 16 to 24 is used to capture the different curvature characteristics of the base plate, side walls, and top arch. Too few control points will lead to underfitting and inability to express local deformation, while too many control points will lead to overfitting and amplification of measurement noise.

[0137] Extract the control point sequence of the fitted B-spline curve as a parameterized contour feature.

[0138] In this embodiment, after the above fitting, the originally hundreds or thousands of random discrete points p l It is compressed and abstracted into a set of ordered control points P j This set of control points constitutes the parametric contour feature. Compared to the original point cloud, the control point sequence has lower data volume and higher signal-to-noise ratio, and possesses clear topological structure information, making it suitable for subsequent deformation transfer analysis.

[0139] In a further embodiment, the method further includes: calculating a residual index for the B-spline fitting, evaluating the fitting quality, and using it as input for subsequent uncertainty analysis. Specifically, after completing the B-spline curve fitting, it is necessary to quantify the degree of deviation between the fitting result and the original measurement data. This deviation reflects the portion of sonar measurement noise that has not been eliminated by the smoothness regularization term. The fitting residual index ε fit Defined as the root mean square value of the distances from all measurement points to the fitted curve. The specific calculation formula is as follows:

[0140] εfit = sqrt((1 / n)*Σ(norm(p l - C(u l )) 2 ));

[0141] Where, ε fit Let be the fitting residual index for this cross-section, n be the total number of measurement points in the cross-section point set, Σ represent the summation of l from 1 to n, and p l Let C(u) be the coordinates of the l-th measured point in the cross-section point set. l ) is a B-spline curve intersecting with p l The corresponding point coordinates, norm represents the Euclidean norm of the vector, and sqrt represents the square root operation.

[0142] Fitting residual index ε fit This reflects the data quality of the cross-section. When ε fit If the value is too large, it indicates that the sonar data for that section is heavily noisy or contains unidentified anomalies, resulting in low reliability of deformation parameters calculated based on that section. The system can set a residual threshold ε. max When ε fit Greater than ε max If this triggers a data quality warning, it is recommended to repeat the scan or perform a manual review of the area.

[0143] In another embodiment of this application, the deformation parameters of the structural seam are calculated based on parametric contour features, including:

[0144] Obtain the sequence of reference control points corresponding to the reference section, and calculate the displacement vector of each measured control point in the control point sequence relative to the reference control point in the corresponding reference control point sequence.

[0145] In other words, the reference control point sequence corresponding to the reference section is obtained from the tunnel reference model. The reference control point sequence adopts the same node vector configuration as the control point sequence. The displacement vector of each measured control point in the control point sequence relative to the corresponding reference control point is calculated.

[0146] Specifically, the system loads the control point sequence P of the reference section. ref The reference cross-section can be an ideal cross-section generated from the tunnel design model, or a B-spline cross-section fitted at the same location during the first or previous inspection. For managing multi-phase inspection data, the system establishes a database indexed by the unique identifier of the structural joint and the inspection phase, storing the control point sequence and its metadata for each phase of inspection. The temporal correspondence between two phases of data is achieved through the following steps: based on the mileage L of the structural joint... seamRetrieve the records of the structural seam in the baseline period from the database; verify whether the B-spline parameter configurations of the two periods are consistent, including the curve order k, the number of control points m+1, and the node vector U. If the configurations are consistent, the control points are directly mapped one-to-one according to the index j; if the configurations are inconsistent, reparameterization is performed before establishing the correspondence. Since a unified node vector configuration was adopted in the aforementioned steps, the measured control point sequence P... meas The j-th point P in meas_j Compared with the j-th point P in the reference sequence ref_j They have a natural correspondence. The system calculates the displacement vector ΔP between each pair of corresponding control points. j ΔP j =P meas_j -P ref_j ; Displacement vector ΔP j It is a three-dimensional vector or a two-dimensional planar vector, comprehensively containing all deformation information of the local region. Through an explicit temporal correspondence mechanism, the calculated displacement vector ΔP j It reflects the real changes of the same physical location at different points in time, rather than spurious displacements caused by differences in parameterization.

[0147] Furthermore, before comparing data from multiple periods, it is necessary to ensure that the B-spline curves fitted from different periods have the same parameterized structure. If the node vector configuration used in the historical benchmark data is inconsistent with the current detection data, one of them needs to be reparameterized to establish a one-to-one correspondence between the two sets of control point sequences. Specifically, by inserting new node values ​​into the original node vectors, the number of control points can be increased without changing the curve geometry. After node insertion, the new control points are calculated from a linear combination of the original control points. For the node u to be inserted... new In the case of u new Located at the original node u j with u (j+1) Between, the new control point P i_new The calculation formula is as follows:

[0148] P i_new = α i * P (i-1) + (1 -α i ) * P i ;

[0149] Among them, P i_new Let P be the i-th new control point after the node is inserted. (i-1) and P i Let α be the original control point. i The mixing coefficient is determined by the relative relationship between the node insertion position and the original node.

[0150] By repeatedly performing node insertion operations, a B-spline curve with arbitrary node vector configuration can be converted into a target node vector configuration, thereby achieving control point alignment across different periods. As a simplified alternative, a uniform node vector configuration can also be enforced across all periods of detection, avoiding the need for reparameterization.

[0151] Obtain the normal vector at each reference control point, and decompose the displacement vector into a radial component along the normal vector direction and a tangential component perpendicular to the normal vector direction.

[0152] Specifically, in engineering and mechanical analysis, displacements in different directions represent different disease mechanisms. Radial displacement typically represents the convergence or expansion of the cross-section (or misalignment for structural joints), while tangential displacement typically represents shear or torsion. The system calculation baseline curve is located at control point P. ref_j The unit normal vector n at the location j For example, the normal vector at the reference control point is obtained by rotating the tangent vector of the B-spline curve. For a planar curve, the normal vector and the tangent vector are perpendicular to each other. The normal vector is obtained by calculating the tangent vector of the curve at the corresponding parameter position of the control point and rotating the tangent vector 90 degrees around an axis perpendicular to the section. Specifically, the tangent vector T(u) of the B-spline curve C(u) at the parameter u is given by the first derivative of the curve with respect to the parameter:

[0153] T(u) = dC(u) / du =Σ(P j * dN j，k (u) / du);

[0154] Where T(u) is the tangent vector at parameter u, dC(u) / du represents the derivative of the curve with respect to parameter u, Σ represents the summation over all control points j, and P j Let dN be the coordinates of the j-th control point. j，k (u) / du is the derivative of the j-th k-th basis function with respect to the parameter u.

[0155] The derivatives of the B-spline basis functions can also be calculated using recursive formulas or approximated by numerical difference methods. After obtaining the tangent vector T(u), the formula for calculating the unit normal vector n(u) is:

[0156] n(u) = R 90 * T(u) / norm(T(u));

[0157] Where n(u) is the unit normal vector at parameter u, R 90 For a two-dimensional plane curve, norm(T(u)) is the rotation matrix that rotates the vector by 90 degrees around the normal axis of the cross section. For a two-dimensional plane curve, norm(T(u)) is the magnitude of the tangent vector.

[0158] Using the principle of vector projection, the displacement vector ΔP j Decompose the components. Radial component ΔP radial_j The calculation formula is:

[0159] ΔP radial_j =(ΔP j *n j )*n j ;

[0160] Where * denotes vector dot product. The magnitude of the radial component represents the amount of stretching or shrinking of the cross-section at that point along the radial direction.

[0161] Tangential component ΔP tangent_j The calculation formula is:

[0162] ΔP tangent_j =ΔP j -ΔP radial_j ;

[0163] The tangential component represents the amount of slippage of the cross section along the circumferential tangential direction.

[0164] The misalignment of the structural joint is calculated based on the radial component, and the opening of the structural joint is calculated based on the tangential component.

[0165] In this embodiment, the misalignment reflects the radial height inconsistency between the pipe sections on both sides of the structural joint. The system selects the difference between the average radial displacement of the section behind the joint and the average radial displacement of the section before the joint as the misalignment d. offset :

[0166] d offset =mean(ΔP after ·n)-mean(ΔP before ·n);

[0167] Where mean represents the average value, ΔP after ·n is the signed projection value of the displacement vector at each control point on the rear section of the structural joint along the direction of the normal vector, ΔP before ·n represents the signed projection of the displacement vector at each control point on the front section of the structural joint along the normal vector direction. Positive values ​​indicate displacement along the normal vector direction (radial outward), and negative values ​​indicate displacement in the opposite direction (radial inward). If the misalignment d offset A positive value indicates that the tube section protrudes after the seam is joined; a negative value indicates that it is recessed.

[0168] The opening reflects the change in distance between the two sides of the structural joint along the tunnel axis. In slice analysis, if three-dimensional control points are considered, the tangential component includes the axial component. Alternatively, the difference in the Z-axis coordinates between the corresponding control points on both sides of the joint can be directly calculated. This is because the normal vector n of each section in the tunnel reference coordinate system... jLocated in a cross-sectional plane perpendicular to the tunnel axis, the component of the displacement vector along the tunnel axis (Z-axis direction) is entirely retained in the tangential component. Therefore, calculating the opening based on the tangential component is equivalent to directly calculating the coordinate difference between the corresponding control points on both sides of the joint in the Z-axis direction:

[0169] d open =mean((z after -z before ) meas -(z after -z before ) design );

[0170] Where d open For opening amount; z after The z-axis coordinate value of the corresponding control point on the rear section of the structural joint is given by z. before Let Z be the Z-axis coordinate of the corresponding control point on the front section of the structural joint. meas Given the measured value, ( ) design To design the ideal value.

[0171] Through decoupled calculations, this embodiment can clearly distinguish whether the structural joint has been misaligned vertically (misalignment), stretched horizontally (opening), or twisted based on the rotation trend of the tangential component, providing accurate diagnostic basis for engineering maintenance.

[0172] This embodiment uses B-spline curves as a mathematical carrier to transform discrete point cloud data into a sequence of control points with physical meaning, and uses vector decomposition technology to achieve qualitative and quantitative analysis of the minute deformation of structural seams.

[0173] In another embodiment of this application, the calculation of the opening amount needs to consider the relative displacement changes of the two phases of detection data in the tunnel axial direction. The opening amount reflects the change in the distance between the pipe sections on both sides of the structural joint in the axial direction; a positive value indicates that the joint has opened up, and a negative value indicates that the joint has compressed.

[0174] Specifically, let the mean axial coordinate of the pipe section before the joint be z. before The mean axial coordinate of the section after the joint is z after The values ​​recorded in the baseline period detection are z. before (1) and z after (1) The values ​​recorded in the current period detection are z. before (2) and z after (2) Opening quantity d open The calculation formula is as follows:

[0175] d open = (z after (2) - z after (1) ) - (z before (2) - z before (1) );

[0176] Where, d open z is the opening amount of the structural joint. after (2) Let z be the mean axial coordinate of the cross-section after the initial suture. after (1) The mean axial coordinate of the cross-section after the reference period is z. before (2) Let z be the mean axial coordinate of the section before the initial seam. before (1) The mean axial coordinates of the cross-section before the reference seam are taken. (1) Indicates base period data, (2) This indicates data from the current period.

[0177] Calculate the axial displacement of the pipe section after the joint and the pipe section before the joint separately. The difference between the two is the change in the width of the structural joint. If the pipe section after the joint moves downstream a greater distance than the pipe section before the joint relative to the reference period, it indicates that the structural joint has been opened.

[0178] In a further embodiment, in addition to misalignment and opening, the structural joint may also undergo relative rotational deformation, i.e., the pipe sections on both sides of the joint undergo relative torsion around the tunnel axis. The relative rotation θ rotation The difference in rotational components of the displacement fields of the control points on both sides of the joint is determined by analyzing the difference in these components. Specifically, for the set of control point displacement vectors {ΔP} in the pipe section region before the joint... j_before The set of control point displacement vectors {ΔP} in the post-joint pipe section region j_after Each rigid body rotation component is fitted. The least squares method is used to solve for the rotation angle that best approximates the displacement field. The rotation component θ of the pre-slit region... before The objective function to be solved is:

[0179] θ before = argmin θ Σ(norm(ΔP j_before - R(θ) * P j_before + P j_before ) 2 );

[0180] Where, θ before argmin is the optimal rotation angle for the displacement field in the pre-suturing region. θLet θ represent the value that minimizes the objective function, Σ represent the summation over all control points j in the pre-suture region, and ΔP represent the summation over all control points j in the pre-suture region. j_before Let P be the displacement vector of the j-th control point in the pre-joint region, R(θ) be the rotation matrix around the tunnel axis by an angle θ, and P be the displacement vector of the j-th control point in the pre-joint region. j_before Let be the coordinates of the j-th control point in the pre-suture region during the baseline period, and norm denote the Euclidean norm of the vector.

[0181] The rotational component θ of the region behind the seam is solved using the same method. after The relative rotation θ of the structural joint rotation The calculation formula is as follows:

[0182] θ rotation = θ after - θ before ;

[0183] The relative rotation θ of the structural joint rotation A positive value indicates that the tube section after the suture rotates clockwise relative to the tube section before the suture, while a negative value indicates counterclockwise rotation.

[0184] In a simplified numerical example, suppose we take one cross-section before and after a structural seam for comparison, and use 8 control points for B-spline curve fitting for each cross-section. Taking control point P0 at the top of the cross-section as an example: [The following is a description of the measured control point P0 on the cross-section before the seam.] 0_before Coordinates [0.000, 2.500, 105.480] meters (X, Y, Z, where Z is the mileage direction); Measured control point P of the post-joint section. 0_after The coordinates are [0.002, 2.508, 105.520] meters; the corresponding reference control point P in the reference model. 0_ref Coordinates are [0.000, 2.500, 105.500] meters; reference control point P 0_ref The unit normal vector n0 at the point is [0.000, 1.000, 0.000], pointing radially outward. Calculate the displacement vector of the control point at the section before the seam: ΔP 0_before = P 0_before - P 0_ref =[0.000, 0.000, -0.020] meters; Calculate the displacement vector of the control point of the section after the joint: ΔP 0_after = P 0_after - P 0_ref =[0.002, 0.008, 0.020] meters; decompose the displacement vector into radial and tangential components. Taking the seam as an example: radial component ΔP radial_after = (ΔP 0_after ·n0)×n0 = (0.008)×[0.000, 1.000, 0.000] = [0.000, 0.008, 0.000] meters; tangential component ΔPtangent_after =ΔP 0_after -ΔP radial_after =[0.002, 0.000, 0.020] meters; calculate the misalignment and opening amount at each control point, and take the average value. Assume the calculation results for the 8 control points are averaged, and the misalignment d... offset = mean(||ΔP radial_after ||) - mean(||ΔP radial_before ||)≈8.0 mm - 0.0 mm = 8.0 mm; opening amount d open ≈0.020 m - (-0.020 m) = 40 mm.

[0185] The results indicate that the structural joint has a radial misalignment of approximately 8 mm and an axial opening of approximately 40 mm. Based on the confidence assessment method, if the calculated confidence interval is ±1.5 mm, the relative uncertainty η = 1.5 / 8.0 ≈ 19%, which is rated as Grade B. It is recommended to make a comprehensive judgment based on historical data.

[0186] According to one aspect of this application, after calculating the deformation parameters of the structural joint, the method further includes a reliability assessment of the detection results of the deformation parameters, specifically including:

[0187] A measurement uncertainty model for sonar ranging is constructed, which describes the functional relationship between single-point ranging error and detection distance and incident angle.

[0188] In this embodiment, the measurement accuracy of the underwater sonar is not a fixed value, but is affected by both the sound wave propagation characteristics and the surface geometry of the object being measured. To accurately quantify the quality of each sampling point, a parameterized error model dependent on the detection distance r and the incident angle θ is established. The measurement uncertainty model can be specifically expressed as the following linear formula:

[0189] σ r =σ0*sqrt(1+α*r+β*tan(θ) 2 );

[0190] Where, σ r σ0 is the standard deviation of single-point ranging, used to characterize the intensity of measurement noise; σ0 is the fundamental error level of the sonar system, which usually depends on the pulse width and sampling frequency of the equipment, for example, a value of 2 mm; r is the radial distance from the sonar transducer to the measuring point on the pipe wall; α is the distance attenuation coefficient, reflecting the decrease in signal-to-noise ratio caused by the divergence of sound waves with distance; θ is the angle between the sound beam and the normal vector of the pipe wall, i.e., the angle of incidence; β is the angle influence coefficient, reflecting the increase in uncertainty caused by the weakening of the echo signal intensity and the elongation of the beam when the angle of incidence increases; tan is the tangent function, and sqrt is the square root function.

[0191] By measuring the uncertainty model, the system assigns an uncertainty label to each sonar point cloud. For example, points on the tube wall directly facing the sonar (θ close to 0) have lower uncertainty, while points on the tube wall edge or structural seam sidewall (θ larger) have increased uncertainty. These differentiated error labels provide the basis for subsequent weighted processing.

[0192] Based on the measurement uncertainty model and the spatial registration process, an error propagation model is constructed to transfer the single-point ranging error to the deformation parameters, and the variance or covariance matrix of the deformation parameters is calculated.

[0193] Specifically, measurement noise propagates layer by layer with the transformation. For example, a three-layer error propagation chain is constructed: a measurement layer, a registration layer, and a parameter layer. In the registration layer, because the input point cloud data contains parameters conforming to σ... r The noise in the distribution causes uncertainty in the calculated spatial transformation matrix (rotation R and translation t'). This uncertainty is expressed using the covariance matrix Σ of the registration parameters. reg To describe.

[0194] As a preferred implementation, Σ is estimated using the inverse of the Hessian matrix at the optimal solution of the registration objective function. reg The Hessian matrix is ​​the matrix of second-order partial derivatives of the registration objective function E with respect to the transformation parameters T, and its mathematical definition is H. ij E equals T i and T j The second-order partial derivatives of the Hessian matrix. The Hessian matrix describes the curvature of the objective function at the optimal solution. When the curvature of the objective function in a certain parameter direction is large, it means that a small deviation of the parameter in that direction will lead to a significant increase in the value of the objective function, thus the estimation of that parameter has high certainty. Conversely, when the curvature is small, the uncertainty of the parameter estimation is large. The inverse matrix of the Hessian matrix is ​​H. -1 The Hessian matrix is ​​proportional to the parameter covariance matrix. In point cloud registration problems, the Hessian matrix can be obtained by the analytical second derivative of the objective function or approximated by numerical difference methods. For point-to-surface ICP objective functions, the Hessian matrix can be expressed as the sum of the outer products of the Jacobian vectors corresponding to all points involved in the registration, and has a closed-form expression.

[0195] For example, since the objective function E is a quadratic approximation of the transformation parameter T, its curvature, i.e., the Hessian matrix H, reflects the stability of the registration result. The larger H is, the sharper the peak and the more determined the parameters. Therefore, the following relationship can be used:

[0196] Σ reg =σ r_avg 2*inv(H);

[0197] Where, σ r_avg The average measurement uncertainty of all points involved in the registration is represented by inv, which denotes the matrix inversion operation.

[0198] In some alternative implementations, if the objective function is too complex to analytically solve for the Hessian matrix, the Monte Carlo method can be used as an alternative. The specific steps are as follows: Add artificially synthesized noise conforming to the measurement uncertainty model to the original point cloud data to generate N sets of perturbed point clouds; perform spatial registration on these N sets of point clouds to obtain N sets of transformation parameter results; statistically analyze the distribution of the N sets of results and calculate their sample covariance matrix as Σ. reg Although the Monte Carlo method is computationally intensive, it is highly versatile and suitable for highly nonlinear scenarios.

[0199] In another possible implementation, an error propagation model is constructed, and the variance or covariance matrix of the deformation parameters is calculated using the following formula:

[0200] Σ def =J f ·Σ reg ·(J f ) T ;

[0201] Where, Σ def Let Σ be the covariance matrix of the deformation parameters. reg J is the covariance matrix of the transformation parameters in spatial registration. f Let J be the Jacobian matrix of the deformation parameters relative to the transformation parameters. f ) T It is the transpose of the Jacobian matrix.

[0202] In this embodiment, deformation parameters such as misalignment and opening are functions of the transformation parameters, denoted as D=f(T). To calculate the covariance matrix Σ of the deformation parameters... def The preferred approach is the first-order Taylor expansion approximation, i.e., the Jacobian transitive method. The covariance matrix Σ... def It is a square matrix, where the diagonal elements correspond to the variances of various deformation parameters (such as misalignment and opening), and the off-diagonal elements correspond to the correlations between parameters. Jacobian matrix J f Σ is the matrix of partial derivatives of the function f(T) with respect to all transformation parameters T in the spatial registration process. Formula Σ def =J f ·Σ reg ·(J f ) TThe final detection accuracy is determined by both registration accuracy and geometric sensitivity. For example, in the Jacobian matrix, if the partial derivative of a deformation parameter with respect to the rotation angle is large, it means that the parameter is very sensitive to rotation errors. Small rotational deviations during registration will be amplified into huge deformation reading errors.

[0203] It should be noted that the Jacobian matrix is ​​the first-order partial derivative matrix of a vector-valued function, used to describe how small changes in the input variables map to changes in the output variables. In error propagation analysis, the Jacobian matrix acts as a linearization amplifier for the propagation of uncertainty from lower to higher levels. Specifically, let the transformation parameter vector D be a function of the transformation parameter vector T, i.e., D = f(T). When T has a small perturbation δT, the corresponding change δD of D can be approximated by a first-order Taylor expansion as follows:

[0204] δD≈J f *δT;

[0205] Where δD is the perturbation vector of the deformation parameter; J f For a Jacobian matrix, its elements (J) f ) ij equals D i For T j The partial derivatives of ; δT is the perturbation vector of the transformation parameters.

[0206] Based on the linear approximation, if the covariance matrix of the transformation parameter T is Σ reg Then the covariance matrix Σ of the deformation parameter D def It can be calculated using the covariance propagation formula. It is the standard result of linear transformation of random variables in probability theory, which is based on the assumption that the perturbation is small enough for the linear approximation to hold.

[0207] The confidence intervals of the deformation parameters are calculated based on the variance or covariance matrix, and the confidence level of the detection results is determined according to the relative proportion of the confidence intervals to the deformation parameter values.

[0208] Specifically, based on the normal distribution assumption, using Σ def The diagonal element is the variance σ D 2 Calculate the standard deviation σ of the deformation parameters D Therefore, a confidence interval (CI) with a specific confidence level can be given:

[0209] CI=[D meas -k'*σ D D meas +k'*σ D ];

[0210] Among them, D measHere, k is the measured value of the deformation parameter; k' is the coverage factor, which is 1.96 for 95% confidence and 2.58 for 99% confidence. For example, the detection result can be expressed as: misalignment = 5.0 mm ± 0.8 mm (95% confidence). This means that the true value has a 95% probability of falling between 4.2 mm and 5.8 mm.

[0211] Based on this, to facilitate quick judgment by non-professionals, a relative uncertainty index η=k'*σ can also be defined. D / abs(D meas This system establishes a reliability grading standard. For example, when η is less than 10%, it is rated as Grade A (high reliability), and the result can be directly used for structural safety assessment or law enforcement. When η is between 10% and 30%, it is rated as Grade B (medium reliability), and it is recommended to make a comprehensive judgment based on historical data or manual review. When η is greater than 30%, it is rated as Grade C (low reliability), indicating that the measurement noise or registration error is too large, masking the true deformation, and it is recommended to conduct encrypted testing again or replace it with high-precision equipment. Through the grading mechanism, engineering risks caused by false alarms or omissions are effectively avoided.

[0212] This embodiment uses an error propagation analysis framework to provide statistically significant scientific basis for engineering maintenance decisions.

[0213] In summary, the underwater detection method for structural joint deformation in pressurized tunnels includes: controlling an underwater robot to travel along the tunnel, adaptively adjusting the travel speed and sonar sampling rate based on real-time joint spacing to acquire high-quality cross-sectional sonar data and positioning data; constructing constraint conditions using geometric features such as the base plate normal and structural joint curvature, performing segmented adaptive spatial registration between measured data and a benchmark model to establish a precise spatial correspondence; parametrically fitting the structural joint region using B-spline curves to extract control point sequence features; calculating the misalignment and opening of the structural joint based on the vector decomposition of control point displacements, and conducting uncertainty assessment.

[0214] In an exemplary embodiment, the process of controlling an underwater robot to move along a tunnel and acquire real-time cross-sectional sonar data and positioning data includes rigorous system initialization calibration and data preprocessing. Specifically, system calibration is a prerequisite for achieving multi-source data fusion. Since sensors such as the cross-sectional sonar, DVL, and inertial navigation unit are physically installed at different locations on the underwater robot, and the coordinate axis directions of each sensor may be inconsistent, rigorous extrinsic parameter calibration is necessary. Let the sonar coordinate system be S, and the underwater robot body coordinate system be B. It is necessary to determine the translation vector T of the sonar center relative to the robot center. sb and rotation matrix R sb During the data acquisition process, the original sonar point P s Point P transformed to volume coordinatesb The calculation formula is: P b =R sb *P s +T sb If this calibration step is ignored, minor installation errors will be amplified into significant geometric distortions during long-distance scanning, such as scanning a straight tube as a spiral tube, causing subsequent spatial registration algorithms to fail.

[0215] After sensor calibration, a dynamic coordinate system transformation is required. Because the underwater robot's attitude sways with the water flow during its movement, each frame of point cloud data acquired by the sonar exists in the instantaneous volume coordinate system at that moment. To integrate point clouds acquired at different times into a unified tunnel coordinate system, dynamic compensation must be performed using real-time attitude information. The tunnel coordinate system is typically defined as a right-handed coordinate system with the tunnel's design axis as the Z-axis and the geoid as the XY plane. Specifically, for any point in the sonar point cloud acquired at time t, the transformation formula from the volume coordinate system to the tunnel coordinate system is as follows:

[0216] p tunnel = R attitude (t) * p body + t position (t);

[0217] Where, p tunnel Let R be the three-dimensional coordinate vector of this point in the tunnel coordinate system. attitude (t) is the attitude rotation matrix constructed from the roll angle φ, pitch angle θ, and yaw angle ψ output by the inertial measurement unit at time t. body Let t be the three-dimensional coordinate vector of this point in the underwater robot's body coordinate system. position (t) represents the position vector of the underwater robot's geometric center in the tunnel coordinate system at time t, obtained by Doppler velocimeter velocity measurement integration or integrated navigation. The attitude rotation matrix R... attitude (t) Following aeronautical convention, the three axes of the volume coordinate system are sequentially rotated by multiplying three basic rotation matrices to align with the tunnel coordinate system around the Z, Y, and X axes. Through a dynamic compensation process, the influence of the underwater robot's attitude swaying on the point cloud geometric accuracy is eliminated, ensuring that the input data for subsequent registration algorithms has a unified spatial reference.

[0218] Furthermore, to address the multipath reflection and reverberation noise unique to underwater sonar data, filtering is performed after acquiring the raw data. A strategy combining statistical filters and radius filters is preferred. The specific implementation process of statistical filtering is as follows: for each point p in the point cloud... i Calculate the average distance d between it and its k nearest neighbors. i Calculate the mean μ of the average distance between all points. dand standard deviation σ d Retain points that meet the following conditions and remove the rest of the outliers: abs(d i -μ d ) <m std *σ d Where abs represents absolute value, m std This is the standard deviation multiple threshold, typically ranging from 1.0 to 2.0.

[0219] Based on statistical filtering, radius filtering is further applied to eliminate false points caused by multipath reflections. The multipath effect refers to the phenomenon where sound waves reflect multiple times between pipe walls before returning to the transducer. Due to the extended propagation path, the system misidentifies these as distant false echo points. Radius filtering checks the number of neighboring points within a specified radius for each point; if the number is too small, it is considered an isolated false point and removed. Specifically, for any point p in the point cloud... i Statistical analysis of its radius r neighbor The number of neighboring points N within the range r (p i By combining statistical filtering and radius filtering conditions, the final set of clean point clouds P is retained. clean The following conditions must be met:

[0220] P clean = { p i ∈P raw :abs(d i - μ d ) <m std *σ d And N r (p i )> n min};

[0221] Among them, P clean P is the filtered clean point cloud set. raw For the original point cloud set, p i Let d be the i-th point in the point cloud. i Let this point and its m std The average distance μ between the nearest neighbors d σ is the population mean of the average distances between all points. d Let m be the standard deviation of the average distance between all points. std The threshold is the multiple of the standard deviation, where abs represents the absolute value, and N r (p i Let p be a point. i At radius r neighbor The number of neighboring points within n min This is the minimum threshold for the number of neighboring points.

[0222] The above joint filtering conditions indicate that an effective measurement point must neither deviate significantly from the population statistically in terms of distance nor have sufficient spatial proximity support. A typical parameter setting is r. neighbor Take two to three times the arc length corresponding to the sonar angular resolution, n min Take 5 to 10 points. This embodiment can effectively filter out isolated noise points caused by suspended impurities, providing clean input data for subsequent B-spline fitting.

[0223] Furthermore, after noise filtering, a measurement uncertainty label needs to be added to each retained measurement point as the data basis for subsequent error propagation analysis. The measurement uncertainty label is a scalar value representing the standard deviation of the ranging accuracy at that point; its magnitude is related to the detection distance and incident angle at that point. Specifically, for any point p in the clean point cloud... i According to its ranging value r i and the angle of incidence θ i Calculate the measurement uncertainty label σ r (p i The calculation formula is as follows:

[0224] σ r (p i =σ0 * sqrt(1 +α* r i +β* tan(θ i ) 2 );

[0225] Where, σ r (p i Let p be a point. i The standard deviation of the ranging; σ0 is the basic ranging accuracy of the sonar system, obtained from the factory calibration of the equipment; r i θ is the radial distance from the transducer to the measuring point; α is the distance attenuation coefficient, reflecting the decrease in signal-to-noise ratio caused by the diffusion of sound wave energy with distance; i θ is the angle between the incident direction of the acoustic beam and the normal vector of the tube wall; β is the angle influence coefficient, reflecting the ranging ambiguity caused by the broadening of the echo when incident at an oblique angle; sqrt represents the square root operation, and tan represents the tangent function.

[0226] Through the pre-labeling process, the quality information of each measurement point is expressed in terms of σ. r (p i The data is explicitly recorded in the form of points, allowing subsequent registration algorithms and deformation calculations to be differentiated and weighted based on the reliability of the points.

[0227] According to one aspect of this application, in actual testing, the system also needs to have the capability to detect and handle abnormal situations. The main abnormal situations and handling strategies include:

[0228] Registration Failure Handling: When iteratively optimizing the registration objective function, if the number of iterations exceeds a preset upper limit and convergence is not achieved, or if the converged objective function value E exceeds a preset failure threshold, the registration of that structural segment is deemed to have failed. In this case, the system automatically marks the data segment as requiring manual review and skips subsequent deformation calculation steps, continuing to process the next structural segment. This segment will be listed separately in the inspection report, and a second scan or manual intervention is recommended.

[0229] Sonar data anomaly handling: During data acquisition, if N consecutive... abnormal If the number of valid points in a frame is below the threshold, it is determined to be a sonar anomaly. The system pauses its movement and attempts to reinitialize the sonar equipment; if the retry fails, an alarm is issued to the operator, and the mileage location where the anomaly occurred is recorded.

[0230] Positioning system anomaly handling: When the DVL (Depth Vault) loses its bottom lock, the system automatically switches to IMU (Immediately Inspector Unit) dead reckoning mode and reduces travel speed to minimize accumulated errors. Simultaneously, the system intensifies monitoring of sonar image features; upon detecting significant features such as structural seams, it immediately performs odometer correction. If the DVL loss time exceeds a preset threshold, the system recommends that operators consider returning to base for re-inspection.

[0231] In another exemplary embodiment, after calculating the deformation parameters of the structural joint, the process further involves transforming abstract numbers into intuitive engineering charts and providing intelligent early warnings based on preset standards. Specifically, to facilitate engineers' comprehensive overview of the health status of long-distance tunnels, tunnel unfolding mapping technology can be employed. This maps the three-dimensional cylindrical tunnel wall onto a two-dimensional plane. The specific mapping logic is as follows: the horizontal axis of the two-dimensional image is defined as the tunnel mileage L, and the vertical axis as the cross-sectional angle θ (0 to 360 degrees). The color value of a pixel corresponds to the local deformation amount, such as radial displacement Δr. The color mapping function can be defined as: Color(L, θ) = Map_Function(Δr(L, θ)); where Map_Function is a function that maps numerical values ​​to a pseudo-color spectrum. For example, blue indicates no deformation, and red indicates an indentation exceeding 20 millimeters. Through the unfolded map, the misalignment of the structural joint will appear as a color abrupt change band across the image, allowing engineers to easily identify high-risk areas for defects.

[0232] Furthermore, the system overlays the benchmark model section (design line) and the measured B-spline fitted section (measured line) onto the same polar coordinate graph or Cartesian coordinate graph, and clearly marks the location of the maximum deformation with arrows or shaded areas.

[0233] The system can also generate trend graphs reflecting the evolution of deformation parameters of the same structural joint over inspection periods, resulting in time-series variation curves. These curves use the inspection date or period number as the horizontal axis and the deformation parameter value as the vertical axis, connecting the points from each inspection to form a line, with error bars representing the confidence interval. Specifically, for a given structural joint, the system retrieves its misalignment sequence {d} from the database across all historical periods. offset (1) d offset (2) , ..., d offset (n)} and the corresponding detection time series {t (1) , t (2) , ..., t (n) Plotting these data points on a time coordinate system allows for a direct observation of whether the deformation parameters remain stable, continue to increase, or exhibit abnormal fluctuations. The time-series variation curve provides a dynamic perspective for structural safety assessment. If the curve shows a stable linear growth trend, extrapolation can be used to predict whether the deformation at a future point in time will exceed the safety threshold; if the curve shows a sudden change or accelerated growth, it may indicate that the tunnel structure is experiencing abnormal conditions, requiring close attention.

[0234] Furthermore, the system incorporates a tiered early warning logic. The triggering condition for an early warning depends not only on the magnitude of the deformation parameter but also on the confidence level. The specific judgment logic is as follows:

[0235] If (d) offset >T danger Furthermore (credibility level == A), a red alert has been triggered, and it is recommended to immediately shut off the water supply for maintenance.

[0236] If (d) offset >T warning Furthermore (credibility level >= B), triggering a yellow alert, it is recommended to include it in the maintenance plan for the following year;

[0237] If (d) offset >T warning However (credibility level == C), a blue warning is triggered, suggesting a second test or manual verification in this area;

[0238] Where, d offset For misalignment, T danger T is the danger threshold. warning This serves as a warning threshold. Through a decision-making mechanism that integrates the magnitude of deformation with data reliability, the false alarm rate and missed alarm rate are minimized.

[0239] This embodiment illustrates the system calibration and noise suppression processing performed at the data acquisition front end to ensure the geometric consistency of the test data and the intuitive usability of the final results, as well as the visualization and early warning mechanisms that should be present at the data output back end. These constitute the auxiliary links of a complete industrial-grade testing system, ensuring an effective closed loop from physical signals to engineering decisions.

[0240] In one detailed embodiment, a water conveyance tunnel has an inner diameter of 3.0 meters and a total length of approximately 2000 meters. It is designed with 200 construction joints spaced approximately 10 meters apart. The current inspection task involves underwater sonar scanning of the section from K100+000 to K102+000 (2000 meters in length), with a focus on the structural joint at K101+055, where a slight misalignment was detected in the previous inspection.

[0241] The system loads a pre-built structural joint mileage distribution database, which records the design mileage of structural joint K101+055 as 1055.0 meters. The underwater robot enters the water at K100+000 to initiate detection. Between K100+000 and K101+050, the robot moves at a speed of v... max Traveling at a cruising speed of 0.5 meters per second, the sonar sampling rate is maintained at f min =2 Hz, cross-sectional sampling interval is approximately 0.25 meters. When the robot moves to K101+050, the system calculates the real-time seam distance d. seam =5.0 meters, entering the deceleration zone. When reaching K101+054.8, i.e., d seam When the distance is approximately 0.2 meters, which is less than the encryption threshold of 0.5 meters, the system begins to reduce its speed and increase the sonar sampling rate. As the underwater robot continues to approach the center of the structural seam, its speed gradually decreases to v. min =0.1 meters per second, the sampling rate is gradually increased to f max =10 Hz. The robot passed through the K101+055 structural seam in a low-speed, high-density mode, acquiring high-density cross-sectional data in the seam area. After passing through the structural seam, the system gradually resumed its cruising speed and continued forward detection.

[0242] After data acquisition, the system performs segmented registration processing on the point cloud data from K101+050 to K101+060. The measured base plate normal vector n is extracted from the base plate region. meas =[0.02, 0.9998, 0.00], and the design normal vector n designThe value [0.00, 1.00, 0.00] has a deviation of approximately 1.1 degrees, which is corrected using the base plate constraint term. The measured structural joint location was identified as mileage 1055.08 meters using curvature abrupt change detection, showing a longitudinal deviation of 0.08 meters from the designed mileage 1055.0 meters, which was corrected using the longitudinal constraint term. The spatial transformation matrix of the pipe segment was obtained after registration convergence. B-spline curve fitting was performed on the cross-sectional point clouds before and after the structural joint to extract the control point sequence. Displacement vector decomposition calculations yielded a misalignment of 6.2 mm and an opening of 12.5 mm for the structural joint. Based on the error propagation model, the 95% confidence interval for this detection result is 6.2 ± 1.0 mm for the misalignment and 12.5 ± 2.0 mm for the opening. The relative uncertainties η are 16% and 16%, respectively, both rated as B (medium confidence). The misalignment of 6.2 mm did not exceed the warning threshold T. warning =10 mm, the system output is normal; however, compared with the previous test record of 4.5 mm misalignment, the deformation has developed, and the system automatically marks it as a concern, suggesting that it be monitored more closely in the next test.

[0243] The inspection results indicate that the K101+055 structural joint is currently under control, but there is a continuous deformation trend, and it is recommended to increase the monitoring frequency. The entire inspection lasted approximately 80 minutes, with about 70 minutes spent on the routine section cruise and about 10 minutes on the intensified structural joint scanning, representing an improvement in efficiency compared to traditional low-speed full-line scanning.

[0244] This invention employs a segmented adaptive registration strategy based on geometric constraints. By dividing a long tunnel into independent structural segments to prevent error accumulation, and using the planar normal vector of the tunnel floor to lock the rotational degrees of freedom around the axis, the longitudinal mileage error is corrected using the curvature abrupt change characteristics of the structural joints. This achieves accurate spatial attitude restoration in an underwater environment lacking external texture. By constructing a prior mileage database and providing real-time distance feedback, the underwater robot is controlled to automatically and smoothly decelerate and increase the sonar sampling rate when approaching the structural joints, and resume cruising speed when moving away from the joint area. This effectively solves the engineering pain points of low efficiency at low speeds and missed features during high-speed scanning. A curve fitting algorithm with a smoothness regularization term is used to effectively suppress the reverberation noise of the underwater sonar. Furthermore, through the normal projection analysis of the control point displacement, the complex mixed deformation is successfully decoupled into radial misalignment and tangential opening, enabling accurate qualitative and quantitative detection of minor defects.

[0245] 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. An underwater method for detecting deformation of structural joints in pressurized tunnels, characterized in that, include: Control the underwater robot to move along the tunnel and acquire real-time cross-sectional sonar data and positioning data; Based on cross-sectional sonar data and positioning data, spatial registration is performed using a pre-stored tunnel reference model to establish a spatial correspondence between the measured cross-section and the reference cross-section. Based on spatial correspondence, parametric contour features of the structural seam region are extracted; The deformation parameters of the structural joint are calculated based on parametric contour features. The deformation parameters include misalignment and opening. Establish the spatial correspondence between the measured cross-section and the reference cross-section, including: Based on the pre-stored structural joint mileage distribution, the cross-sectional sonar data is divided into predetermined independent structural segment point clouds; For each structural segment point cloud, a registration objective function is constructed, which includes a data fitting term to minimize the point cloud matching error and a geometric feature constraint term to limit the degrees of freedom of rigid body transformation. By iteratively optimizing the registration objective function, the spatial transformation matrix that minimizes the registration objective function is obtained, which serves as the spatial correspondence. Extract parametric contour features of the structural seam region, including: By utilizing spatial correspondence, the cross-sectional sonar data is unified to the tunnel reference coordinate system, and the point cloud of the structural joint area is sliced ​​along the tunnel axis to obtain a predetermined set of cross-sectional points. For each cross-sectional point set, B-spline curve fitting is performed, and a pre-defined uniform node vector configuration is used during the fitting process; Extract the control point sequence of the fitted B-spline curve as a parameterized contour feature; Deformation parameters of structural joints are calculated based on parametric contour features, including: Obtain the sequence of reference control points corresponding to the reference section, and calculate the displacement vector of each measured control point in the control point sequence relative to the reference control point in the corresponding reference control point sequence. Obtain the normal vector at each reference control point, and decompose the displacement vector into a radial component along the direction of the normal vector and a tangential component perpendicular to the direction of the normal vector. The misalignment of the structural joint is calculated based on the radial component, and the opening of the structural joint is calculated based on the tangential component.

2. The method according to claim 1, characterized in that, After acquiring real-time positioning data, the process also includes dynamically adjusting the underwater robot's travel speed and sonar sampling rate, specifically: Based on positioning data, the real-time gap between the underwater robot’s current position and the nearest structural gap in the pre-stored list of structural gap positions is calculated. The underwater robot's travel speed and sonar sampling rate are dynamically adjusted based on the real-time gap. When the real-time gap is less than the preset encryption threshold, the travel speed is reduced and the sonar sampling rate is increased.

3. The method according to claim 2, characterized in that, The underwater robot's travel speed and sonar sampling rate are dynamically adjusted based on the real-time gap. Specifically, the target travel speed and target sampling rate of the underwater robot are calculated using the following formula: v(t)=v max -(in max -v min )*exp(-(d seam (t)) 2 / (2*σ v 2 )); f s (t)=f min +(f max -f min )*exp(-(d seam (t)) 2 / (2*σ f 2 )); Where v(t) is the target speed, f s (t) represents the target sampling rate, d seam (t) represents the real-time seam distance, v max and v min These are the preset maximum and minimum speeds, f. max and f min These are the preset maximum and minimum sampling rates, σ v and σ f The preset adjustment range parameter is exp(), which is an exponential function.

4. The method according to claim 1, characterized in that, Geometric feature constraints include base plate plane constraints, and their construction process includes: Extract the base plate region point set from the structural segment point cloud, and use a plane fitting algorithm to calculate the measured base plate normal vector of the base plate region point set; Obtain the corresponding design floor normal vector from the tunnel reference model; The constraints on the bottom plate plane are constructed to minimize the angle between the measured bottom plate normal vector and the designed bottom plate normal vector, and the rotational degrees of freedom about the tunnel axis are locked.

5. The method according to claim 4, characterized in that, The geometric feature constraints also include longitudinal position correction constraints, the construction process of which includes: Calculate the rate of change of profile curvature of each section in the point cloud of the structural segment, and identify the location where the rate of change of profile curvature exceeds the preset abrupt change threshold as the actual structural seam location; Obtain the design structural joint position corresponding to the measured structural joint position, and calculate the longitudinal deviation between the two. Constraints on the longitudinal position are constructed based on the longitudinal deviation value, the mileage component in the positioning data is corrected, and the translational degree of freedom along the tunnel axis is locked.

6. The method according to claim 5, characterized in that, The registration objective function is solved through iterative optimization. Specifically, the minimum value of the following registration objective function is found to determine the rotation matrix R and the translation vector t': E=Σ||R·p i +t'-q i || 2 +λ1·||R·n meas -n design || 2 +λ2·||t z -ΔL|| 2 ; Where E is the registration objective function value, p i For points in the structural segment point cloud, q i For the tunnel reference model and p i The corresponding nearest point, n meas Let n be the measured base plate normal vector. design To design the base plate normal vector, t z Let t' be the component of the translation vector along the tunnel axis, ΔL be the longitudinal deviation value, λ1 and λ2 be the preset constraint weight coefficients, and ||| represent the Euclidean norm.

7. The method according to claim 1, characterized in that, After calculating the deformation parameters of the structural joint, the reliability assessment of the detection results of the deformation parameters is also included, specifically: A measurement uncertainty model for sonar ranging is constructed to describe the functional relationship between single-point ranging error and detection distance and incident angle; Based on the measurement uncertainty model and the spatial registration process, an error propagation model is constructed to transfer the single-point ranging error to the deformation parameter and calculate the variance or covariance matrix of the deformation parameter. The confidence intervals of the deformation parameters are calculated based on the variance or covariance matrix, and the confidence level of the detection results is determined according to the relative proportion of the confidence intervals to the deformation parameter values.

Images