A mine shaft wall health detection method based on multi-source data fusion

Abstract

A kind of vertical shaft wall health detection method based on multi-source data fusion, through time alignment, constructing shaft wall geometry prior factor, constructing cross-modal consistency factor, constructing IMU pre-integration factor, constructing external parameter slow change prior factor, constructing odometry factor, sliding window joint optimization, point cloud data multi-modal disease parallel detection, fusion structured light, laser and IMU multi-source information, realize the vertical shaft wall high-precision three-dimensional reconstruction and disease detection.The present application can still keep the stable space relationship of multi-sensor and the consistency of point cloud splicing under complex working conditions such as vibration, water film reflection, dust shielding, etc., while ensuring the stability, precision and robustness, significantly inhibiting the external parameter drift and point cloud distortion, solving the measurement discontinuity and deformation evaluation distortion problem caused by data quality fluctuation in traditional single or relative splicing scheme, providing technical support for accurate acquisition of shaft wall shape and disease detection.

Description

Technical Field

[0001] This invention relates to a method for detecting the health of vertical shaft walls based on multi-source data fusion, belonging to the technical fields of three-dimensional measurement, sensor fusion positioning and underground engineering structure health monitoring. Background Technology

[0002] As a critical infrastructure in mining and underground engineering, the structural safety of vertical shafts directly impacts mine production safety and personnel safety. During long-term operation, the shaft wall structure may be affected by various factors such as changes in ground pressure, surrounding rock deformation, groundwater seepage, temperature variations, and construction disturbances, resulting in deformation phenomena such as cracking, bulging, convergence, and misalignment. Therefore, conducting high-precision, long-term, and stable deformation monitoring of vertical shaft walls is of great significance for ensuring the safe operation of mines.

[0003] Existing wellbore monitoring methods largely rely on manual measurement, contact sensors, or single-sensor 3D scanning technology. For example, traditional total station measurements and manual inspections are inefficient and difficult to implement continuous high-density monitoring. While non-contact methods based on single laser scanning or visual measurement can acquire wellbore point cloud data, they are susceptible to factors such as dust, water mist, reflective water films, and insufficient lighting in the downhole environment, leading to data loss or increased noise, thus affecting measurement accuracy and stability. Furthermore, existing 3D reconstruction methods often rely on relative point cloud stitching and registration, which can easily generate cumulative errors during long-distance vertical shaft scanning, causing overall model drift and making it difficult to achieve absolute quantitative assessment of wellbore deformation. Simply relying on inertial navigation or odometers also suffers from long-term drift problems and cannot provide a high-precision, stable spatial reference on its own.

[0004] With the development of smart mines and digital twin mines, higher requirements are placed on the accuracy, data stability, and automation level of shaft wall morphology acquisition. Therefore, there is an urgent need for a shaft wall health monitoring method that can operate stably in complex underground environments while balancing accuracy and robustness. This method should improve the reliability and traceability of shaft wall deformation monitoring through the complementarity and fusion of information from multiple sensors. Summary of the Invention

[0005] The purpose of this invention is to provide a method for detecting the health of vertical shaft walls based on multi-source data fusion. This method can work stably in complex downhole environments, balancing accuracy and robustness, improving the reliability and traceability of well wall deformation monitoring, and providing technical support for accurate acquisition of well wall morphology and detection of defects.

[0006] To achieve the above objectives, the present invention provides a method for detecting the health of a vertical shaft wall based on multi-source data fusion, comprising the following steps:

[0007] S1. Time alignment: Establish a unified time reference, map the timestamps of the laser scanner, structured light camera and odometer to the unified time reference, and interpolate the mapped observations to obtain synchronous observations on the same time axis as the IMU status.

[0008] S2. Constructing prior geometric factors for the wellbore: Based on the spatial distribution of laser point clouds in the wellbore reference coordinate system, fit the geometric model of the cylindrical surface of the wellbore and construct the geometric prior residuals from the points to the geometric model;

[0009] S3. Construct cross-modal consistency factor: After transforming the laser point cloud and structured light point cloud to a unified coordinate system, slice or layer them according to the axial depth of the well shaft, and construct cross-modal consistency residuals for the cross-sectional parameters of the two modal point clouds within the same depth slice.

[0010] S4. Constructing IMU pre-integration factors: Based on the angular velocity and acceleration measurements of the IMU between adjacent keyframes, pre-integration is performed to construct the kinematic residuals between the states of adjacent keyframes.

[0011] S5. Constructing Slowly Varying Extrinsic Parameter Prior Factors: Set slowly varying prior residuals for the extrinsic parameters between the laser scanner and the IMU, and between the structured light camera and the IMU, to constrain the variation of the extrinsic parameters within the sliding window;

[0012] S6. Constructing the odometer factor: Based on the absolute rope length output by the odometer, construct the axial displacement residual to provide a dimensional constraint on the axial depth of the wellbore.

[0013] S7. Sliding window joint optimization: The pose, velocity, IMU zero bias, extrinsic parameters, and time offset within the window are the variables to be estimated. The residuals from steps S2 to S6 and the marginalized prior residuals are combined to form the objective function and iteratively minimized to obtain the updated pose, extrinsic parameters, and time offset.

[0014] S8. Parallel detection of multimodal defects using point cloud data: Utilizing the high density and sensitivity to fine surface undulations of structured light point clouds for crack detection; leveraging the geometric stability of laser point clouds and their minimal impact from surface reflection and localized pores, combined with the abnormal responses of structured light point clouds, such as increased porosity under humid conditions, to determine humid or leaking areas.

[0015] Furthermore, the specific process of S1 is as follows:

[0016] S1.1 Unified Clock: A unified clock source is used for timestamp calibration of the laser scanner, structured light camera, inertial measurement unit, and odometer. The unified clock source is the same master control clock distribution + hardware trigger input, ensuring that the data of the laser scanner, structured light camera, IMU and odometer have the same time reference.

[0017] S1.2 Software Time Alignment: Based on a unified clock, using the IMU as a high-frequency reference, the observation timestamps of the structured light camera, laser scanner, and odometry are aligned, mapping the observations to the time intervals between two adjacent IMU frames. The IMU states are then interpolated to obtain the pose at the corresponding time for subsequent residual construction. The process is described as follows:

[0018] ;

[0019] ;

[0020] ;

[0021] in, These are the timestamps from the structured light camera, laser scanner, and odometer, respectively. These represent the time deviations of the observation times of the structured light camera, laser scanner, and odometer relative to the IMU clock. This indicates the corresponding time point on the IMU timeline for observations by the structured light camera. This indicates the corresponding moment on the IMU timeline observed by the laser scanner. This indicates the corresponding moment on the IMU time axis for the odometer observation.

[0022] Furthermore, the specific process of S2 is as follows:

[0023] S2.1 Constructing the cylindrical surface model:

[0024] S2.1-1, in the window Within, for each keyframe Take the corresponding laser point Then, use the pose and extrinsic parameters of that keyframe to transform the laser point into the world frame:

[0025] ;

[0026] in, Let the laser point be located in world coordinates at the k-th keyframe within the window. for The pose transformation from the body coordinate system to the world coordinate system at any given moment. This is an extrinsic parameter transformation from the laser coordinate system to the body coordinate system. Let be the laser point in the laser coordinate system at the k-th keyframe;

[0027] The resulting window point set is obtained by aggregation:

[0028] ;

[0029] in, This refers to all laser points within window W in the world coordinate system.

[0030] S2.1-2. Layer according to the z-axis to obtain a set of slices. :

[0031] ;

[0032] in, Number the slices. Let z be the z-coordinate of the laser point in the world coordinate system within window W. The minimum z-coordinate value of the laser point in the world coordinate system within window W. The thickness of the slice. This is the m-th slice;

[0033] S2.1-3. In each slice, project the point cloud onto a plane perpendicular to the axis and fit a circle using the RANSAC algorithm: output the center point of each slice. , The average z-coordinate of the slice is given by the radius of the circle. Perform axis fitting and set the center point of each layer. Perform linear fitting:

[0034] ;

[0035] The axis to be fitted is represented as:

[0036] ;

[0037] in, For a point on the fitted axis, It is a direction vector. The scale value;

[0038] point The square of the shortest distance to the line is:

[0039] ;

[0040] Perform least-squares line fitting:

[0041] ;

[0042] The average of the fitted radii for each layer is used as the radius of the window. Finally, the cylindrical parameters of the window are obtained. ;

[0043] S2.2 Introduction of slice quality gating: Given that downhole dust, water film reflection, and occlusion factors cause sparse laser point clouds or the presence of anomalous echoes, to avoid unreliable fitting misleading joint optimization, slice quality gating is introduced for cylindrical fitting, forming a slice set participating in geometric priors. ; Calculate the fit quality index for each slice:

[0044] S2.2-1, Interior Point Ratio:

[0045] ;

[0046] in, This represents the number of interior points when fitting a circle using the RANSAC algorithm. This represents the number of point clouds within the slice; , Take a value of 0.5 to 0.8;

[0047] S2.2-2, Coverage:

[0048] Project the slice points onto the cross-sectional plane, calculate the polar angle distribution of the points, and obtain the coverage angle. :

[0049] ;

[0050] in, , Take a value of 0.3 to 0.5;

[0051] S2.2-3, Fitting Residual: The root mean square error of the fit for all points is:

[0052] ;

[0053] ;

[0054] in, For this slice of point cloud set, For a point in this slice set, Let be a point on the axis. This represents the radial distance residual from the point to the cylindrical surface; , This is laser ranging noise (as specified by the manufacturer). Values ​​range from 2 to 5;

[0055] S2.2-4, when satisfied When, the slice is included If the geometric prior factors are not constructed, then the geometric prior factors for that slice are constructed.

[0056] S2.3 Computational geometric prior residuals: cylinder parameters For any point The axis is formed by Definition, point Shortest distance to the axis:

[0057] ;

[0058] ;

[0059] in, , The parameters of the cylinder model obtained by fitting within the window include the cylinder axis. and cylinder radius , Point Relative to cylindrical model parameters The geometric prior residuals;

[0060] S2.4, For each Take the point cloud set within the slice. Constructing geometric prior residuals:

[0061] ;

[0062] in, For the corresponding Point cloud within a slice.

[0063] Furthermore, the specific process of S3 is as follows:

[0064] S3.1 Calculate the surface-to-surface residual vector: within the depth range covered by the corresponding window. First, the laser point cloud and structured light point cloud are preprocessed with noise reduction, filtering, and coordinate unification. Then, the point cloud is divided into several depth slices according to the shaft axial depth, forming a depth slice set.

[0065] ;

[0066] in, Let j represent the j-th depth slice, and n represent the number of depth slices.

[0067] For slice unit Point set within The parameters are as follows:

[0068] S3.1-1 Calculate the centroid:

[0069] ;

[0070] Where N is the number of point clouds involved in the calculation in the corresponding slice;

[0071] S3.1-2, Calculate the covariance:

[0072] ;

[0073] S3.1-3, Perform eigenvalue decomposition :

[0074] Eigenvector corresponding to the smallest eigenvalue That is, the plane normal:

[0075] ;

[0076] in, Let be the i-th eigenvalue of the covariance matrix S. For eigenvalues The corresponding feature vector;

[0077] S3.1-4, Planar Offset:

[0078] ;

[0079] Plane offset difference directly reflects the distance difference along the normal direction. It is not enough for the normal directions to be consistent; the two planes may be parallel but offset. Used to constrain translation difference along the normal direction;

[0080] S3.1-5, Normal Consistency Residual :

[0081] ;

[0082] in, Let be the unit normal vector of the locally fitted plane of the laser point cloud. Let be the unit normal vector of the locally fitted plane of the structured light point cloud. This represents the cross product operation of vectors; if the two normals are parallel (in the same direction or in opposite directions), the cross product result is 0, i.e. It is a three-dimensional vector residual;

[0083] S3.1-6, Consistent planar position (Directly reflects the distance difference along the normal direction):

[0084] The same plane can be written as or In comparison First, ensure that the two normal directions are consistent; otherwise... The inversion of the sign will cause the residual to jump. and Let be the plane offsets of the local fitting planes of the laser point cloud and the structured light point cloud, respectively, and let be the plane offsets of the local fitting planes of the laser point cloud and the structured light point cloud.

[0085] ;

[0086] If not satisfied, then:

[0087] ;

[0088] ;

[0089] S3.1-7, Surface-to-Surface Residual Vector:

[0090] ;

[0091] in, This represents the cross-modal surface-to-surface consistency residual vector between the laser point cloud and the structured light point cloud within the corresponding slice z.

[0092] S3.2, Apply confidence weight to the residuals: For extrinsic parameters sensitive to rotation or translation, if a mode mismeasures a surface in a certain locality, it will generate a large erroneous gradient, directly updating the extrinsic parameters in the incorrect direction. To address the issue of residual vectors automatically weakening when data differences occur, this invention calculates confidence weights based on the cross-sectional quality of structured light point clouds and laser point clouds within each depth slice. This allows slices with poor data quality to automatically reduce their contribution during the optimization process.

[0093] ;

[0094] in, Let the cross-membrane state confidence weight be the value at depth j. The cross-sectional confidence level of the structured light point cloud within the j-th slice; Let the cross-sectional confidence level of the laser point cloud within the j-th slice be denoted as .

[0095] Based on the aforementioned normal consistency residual Consistency residuals with planar position The combination yields the first Cross-modal consistency residual vectors for depth slices:

[0096] ; in, Indicates the first The cross-modal consistency residual vector between the laser point cloud and the structured light point cloud within a depth slice; This represents the normal consistency residual of the local fitting plane of the two modal point clouds within the slice; This represents the planar position consistency residual of the local fitting plane of the two modal point clouds within the slice.

[0097] Obtaining transmembrane consistent residuals Then, it is weighted and represented as:

[0098] ;

[0099] in, The transmembrane consistency weighted residual corresponding to the j-th depth slice; Let be the transmembrane consistency residual vector between the structured light point cloud and the laser point cloud in the j-th depth slice; For robust functions, used to reduce the impact of abnormal residuals on the results;

[0100] The credibility formula is:

[0101] ;

[0102] For each slice Projecting the point onto the cross-sectional plane (and) (Orthogonal), fitting circle parameters to structured light slices Fitting circle parameters to laser slices The fitting method uses RANSAC + least squares, where:

[0103] (1) The in-point ratio when fitting the cross section in RANSAC:

[0104] ; Let be the number of interior points. This represents the total number of point clouds;

[0105] use The function avoids the influence of extreme values ​​and maps to [0,1]:

[0106] ;

[0107] (I basically don't believe anything below this).

[0108] (Above this is considered good enough);

[0109] in, This indicates the quality index to be normalized. This indicates the lower confidence limit of the quality indicator, when... When the normalization result is zero, the result is taken as 0. This indicates the upper limit of confidence for the quality indicator, when When the normalization result is 1, the result is taken as 1. This means that the calculation result will be limited to the range [0,1].

[0110] (2) The root mean square distance from the point to the fitted circle after cross-section fitting is:

[0111] ; It is the radial error from the point to the cross section;

[0112] (3) Parameters This is a scaling parameter (determining the descent rate), taken as a multiple of the sensor's ranging standard deviation, to reflect the rapid decay of reliability when the fitting residual exceeds the noise level.

[0113] Laser point cloud: ;

[0114] Institutional Light: ;

[0115] in, This represents the standard deviation of the ranging noise of the laser scanner. The standard deviation of the ranging noise of the structured light camera;

[0116] (4) The coverage angle is determined by whether the slice is sufficiently covered (to prevent fitting a circle even when only a small arc is visible). This indicates the degree of coverage of the slice on the cross-section, i.e., the angular coverage rate. By projecting the slice points onto the cross-sectional plane and calculating the polar angle distribution of the points, the coverage angle can be obtained. :

[0117] ;

[0118] use To avoid the influence of extreme values, map to [0,1]:

[0119] .

[0120] (Approximately 72°C; below this temperature, degradation is likely).

[0121] (Approximately 216°, above which it is basically stable).

[0122] in, This represents the angle coverage to be normalized. This represents the lower bound of the coverage confidence level, when At that time, it was considered that the slice coverage was insufficient, and the normalization result was set to 0; Indicates the upper limit of coverage confidence, when At that time, the slice was considered to have sufficient coverage, and the normalization result was set to 1; This means that the value is restricted to the range [0,1].

[0123] Furthermore, the specific process of S4 is as follows:

[0124] S4.1 Determine the parameters to be estimated;

[0125] S4.1-1, Window Status:

[0126] ;

[0127] in, These represent IMU's position in the world system. The following posture, position, and speed; These represent zero bias of the IMU gyroscope and zero bias of the accelerometer, respectively.

[0128] S4.1-2, Global parameters to be estimated:

[0129] ;

[0130] in, These represent the extrinsic parameters from the IMU to the camera and the extrinsic parameters from the IMU to the laser (rigid body transformation), respectively. These represent the offsets of the camera, laser, and odometry timestamps relative to the IMU, respectively.

[0131] Connecting adjacent states within the window provides constraints, while zero bias is introduced into the pre-integral model for estimation. ;

[0132] The IMU pre-integration residuals include rotational residuals, velocity residuals, position residuals, gyro zero-bias residuals, and accelerometer zero-bias residuals.

[0133] ;

[0134] S4.2 Parameter Definition: Select a set of keyframes within the sliding window. For each adjacent keyframe An IMU pre-integration factor is constructed for each frame to constrain the relative motion between adjacent keyframes:

[0135] S4.2-1, Set adjacent keyframes The corresponding timestamps are respectively Since the IMU sampling frequency is usually higher than the keyframe frequency, there are multiple IMU sampling times between two adjacent keyframes, denoted as:

[0136] ;

[0137] in, Keyframe IMU sampling times between; This is expressed as the number of IMU samples within this adjacent keyframe interval, and satisfies:

[0138] ;

[0139] If the keyframe timestamp does not coincide with the IMU sampling time, time interpolation is used to obtain the IMU status or measurement value at the endpoint;

[0140] Initialize the pre-integral quantity.

[0141] ;

[0142] in, For initialization from keyframes arrive Pre-integral relative rotation; Represents a third-order identity matrix; This is the pre-integral relative velocity increment during initialization; This represents the pre-integration relative position increment during initialization; 0 indicates the zero vector.

[0143] S4.2-2, Keyframes arrive At the m-th IMU sampling time, the gyroscope angular velocity measurement and accelerometer measurement values ​​output by the IMU are denoted as follows: and Remove the zero bias from it:

[0144] ;

[0145] in, Indicates the first The angular velocity vector after removing the gyroscope zero bias at the m-th IMU sampling time within a keyframe interval; This represents the acceleration vector after removing the accelerometer zero bias at the m-th IMU sampling time within the i-th keyframe interval; This indicates that the gyroscope has zero bias. This indicates that the accelerometer has zero bias.

[0146] The vector representation of the angular velocity after zero bias is:

[0147] ;

[0148] in, , , These represent the deflection angular velocity at the m-th IMU sampling time within the i-th keyframe interval. , , Components in three directions;

[0149] S4.2-3, Rotational pre-integral update:

[0150] For each IMU sampling interval Perform rotational pre-integral update:

[0151] ;

[0152] The time interval between two adjacent IMU sampling times is defined as:

[0153] ;

[0154] in, Indicates from keyframe The cumulative pre-integral relative rotation up to the m-th IMU sampling time; express The exponential mapping on the top is used to convert rotation vectors into rotation matrices; during initialization... As the IMU sampling data is gradually integrated, This indicates that the cumulative relative rotation is no longer equivalent to the identity matrix;

[0155] S4.2-4, Velocity pre-integral update is as follows:

[0156] ;

[0157] in, To start from keyframe The cumulative relative velocity increment obtained up to the m-th IMU sampling time; at this time Used to transfer a vector from the current time-space system Switch to the initial frame machine system ;

[0158] S4.2-5, Position pre-integration updated as follows:

[0159] ;

[0160] in, Indicates from keyframe The pre-integrated relative position increment accumulated up to the m-th IMU sampling time; This represents the debiased acceleration under the sampling time of the m-th IMU. Transition to the starting keyframe Under the machine system, thus accumulating velocity and position pre-integral quantities in the same coordinate system;

[0161] Complete keyframes and Integrating over all IMU sampling intervals, we obtain the IMU pre-integration between adjacent keyframes:

[0162] ;

[0163] ;

[0164] ;

[0165] in, , and These represent keyframes. To keyframe The relative rotation, relative velocity increments, and relative position increments are obtained by pre-integration of the IMU sampled data.

[0166] S4.3 Residual calculation in:

[0167] S4.3-1 Rotational residuals in:

[0168] ;

[0169] , ;

[0170] It is a logarithmic mapping of rotation matrices, outputting a 3D rotation vector;

[0171] This represents the relative rotation prediction obtained by integrating the gyroscope measurements after zero bias removal during the period from i to i+1. This represents the actual rotation obtained by integrating the gyroscope measurements after removing the zero bias during the period from i to i+1;

[0172] S4.3-2 Velocity residual in the middle:

[0173] ;

[0174] in, and These are represented as keyframes. and The velocity of the IMU in the world coordinate system at any given moment; Represented as the gravitational acceleration vector in the world coordinate system; Represented as keyframe and The time interval between; This indicates that the keyframe is estimated from the current state and then transformed. The relative velocity increment under the machine system; This represents the relative velocity increment obtained from IMU pre-integration;

[0175] S4.3-3 Positional residuals in:

[0176] ;

[0177] in, and These are represented as keyframes. and The position of the IMU in the world coordinate system at any given moment; Keyframe The velocity of the IMU in the world coordinate system at any given moment; This is represented as an estimate obtained from the current state and transformed into a keyframe. The relative position increment under the machine system; Represented as the relative position increment obtained from IMU pre-integration;

[0178] S4.3-4 Zero-biased residuals in:

[0179] , ;

[0180] in, The gyroscope has zero bias residual; The accelerometer has zero bias residual; and These are represented as keyframes. and The gyroscope is at zero bias at any given moment; and These are represented as keyframes. and The accelerometer zero bias at time 1; this zero bias residual is used to constrain the continuity of the IMU zero bias between adjacent keyframes, avoiding unreasonable jumps in the zero bias estimation;

[0181] S4.4, IMU residual vector:

[0182] Construct adjacent keyframe pairs based on keyframes within the sliding window. :

[0183] ;

[0184] in, N-1 is the set of adjacent keyframe pairs that participate in the IMU pre-integration constraint within the sliding window; N-1 is the number of keyframes within the sliding window.

[0185] For any pair of adjacent keyframes ,in By combining the rotational residual, velocity residual, position residual, gyroscope zero-bias residual, and accelerometer zero-bias residual, the IMU pre-integration residual vector is obtained:

[0186] .

[0187] Furthermore, the specific process of S5 is as follows:

[0188] S5.1, Marginalized Prior Residuals: Marginalized Prior Residuals As the sliding window moves forward, the oldest state variable that moves out of the window is marginalized. Its constraint on the objective function is compressed into a priori factors and applied to the retained variables within the window, forming marginalized prior residuals, expressed as:

[0189] ;

[0190] in, To preserve the increment of variables within the window, the matrix with vector It is obtained by marginalization elimination (Schur complement), which makes the cost function before and after marginalization equivalent in the reserved variable space;

[0191] S5.2, Slowly Varying Prior Residue of Extrinsic Parameters: To suppress abrupt changes in extrinsic parameters caused by misconstraints due to downhole vibration, temperature drift, and point cloud outliers, slowly varying prior residuals are introduced into the extrinsic parameters. The keyframe set is maintained within a sliding window. The body to laser external parameters Define the set of adjacent keyframe pairs within the window:

[0192] ;

[0193] For any adjacent keyframe pair Define the extrinsic parameter as a slowly varying prior residual:

[0194] ;

[0195] in, The rotation matrix represents the attitude of the laser coordinate system relative to the body coordinate system. The translation vector represents the position of the origin of the laser coordinate system within the body coordinate system.

[0196] Exterior parameters from the body to the structured light camera ,definition:

[0197] ;

[0198] in, for Logarithmic mapping is used to represent rotational differences as three-dimensional vectors;

[0199] The extrinsic parameter estimates at time i are combined to obtain the saturated prior of the extrinsic parameters:

[0200] ;

[0201] in, These are the covariance matrices for the changes in extrinsic parameters of laser and structured light, used to allow for slow drift of extrinsic parameters and suppress non-physical jumps; preferably, the standard deviation of the rotation increment is taken as... The standard deviation of the translation increment is taken as ;

[0202] S5.3, Prior Residual Term It is defined as a combination of marginalized prior residuals and slowly varying extrinsic prior residuals.

[0203] ;

[0204] in, To balance the relative strengths of the two types of priors, the extrinsic slowly varying prior weight coefficients are used. Reliability of cross-modal constraints Related settings; credibility Credibility of cross-modal sections of slices at various depths within the sliding window In summary, the slowly varying prior weights of the extrinsic parameters are adaptively adjusted according to the following relationship:

[0205] ;

[0206] in, As the benchmark weight, The adjustment coefficient controls the magnitude of the increase in prior weights when confidence decreases. Both can be set or adjusted online based on equipment rigidity, downhole vibration and temperature drift levels, and data quality indicators.

[0207] Furthermore, the specific process of S6 is as follows:

[0208] S6.1 Calculate the odometer increment:

[0209] The odometer outputs the absolute distance for any IMU at critical moments. The odometer increment is defined as:

[0210] ;

[0211] in, This represents the axial displacement increment obtained from the odometer between keyframe i and keyframe j. This represents the absolute distance reading of the odometer at time t; It represents the time offset of the odometer timestamp relative to the IMU timestamp, and is used to map critical moments of the IMU to the odometer time axis;

[0212] Because the odometer sampling frequency may be lower than that of the IMU, when the required time When the odometer sampling time does not coincide with the odometer sampling time, linear interpolation is used to obtain the odometer reading at that time. ,but:

[0213] ;

[0214] in, and The sampling times of two adjacent odometers; Indicates the odometer at time The absolute distance reading, Indicates the odometer at time The absolute distance reading, Indicates the time obtained through linear interpolation. The corresponding absolute distance reading from the odometer;

[0215] S6.2 Calculate the axial displacement increment estimated by the IMU. :

[0216] Based on the IMU pose obtained through sliding window optimization, calculate the displacement increment along the well shaft axis between two IMU keyframe moments:

[0217] set up and They represent and The position vector of the IMU system origin in the world coordinate system at any given time. Let represent a predefined unit vector along the wellbore axis. Then, the axial displacement increment estimated by the IMU is:

[0218] ;

[0219] in, This represents the unique increment of the wellbore axis obtained by the pose estimation of the IMU between keyframe i and keyframe j; They are respectively Displacement estimated by the IMU at time step; The defined unit vector along the wellbore axis; This indicates that the three-dimensional displacement between two keyframes is projected onto the shaft axial direction.

[0220] S6.3 Calculate the odometer residual:

[0221] Define the set of adjacent keyframe pairs within the window:

[0222] ;

[0223] For any adjacent keyframe pair Define the extrinsic parameter as a slowly varying prior residual:

[0224] .

[0225] Furthermore, the specific process of S7 is as follows:

[0226] S7.1, Parameters to be estimated:

[0227] Without a reference loop, extrinsic parameters are difficult to obtain stably using a single observation. Within a sliding window, IMU pose, velocity, position, IMU bias, sensor extrinsic parameters, and the time offset of each sensor relative to the IMU are used as joint parameters to be estimated. Through IMU motion constraints, wellbore geometry priors, cross-modal consistency constraints, slowly varying extrinsic parameter priors, and odometer residual constraints, the extrinsic parameters are continuously pulled back to correct values ​​while avoiding miscalibration. This ensures that IMU pose optimization and extrinsic parameter updates are performed within the same optimization framework. The parameters to be estimated consist of two parts:

[0228] S7.1-1, Window Status:

[0229] ;

[0230] in, These represent IMU's position in the world system. The following posture, position, and speed; These represent zero bias of the IMU gyroscope and zero bias of the accelerometer, respectively.

[0231] S7.1-2, Global parameters to be estimated:

[0232] ;

[0233] in, These represent the extrinsic parameters from the IMU to the camera and the extrinsic parameters from the IMU to the laser (rigid body transformation), respectively. These represent the offsets of the camera, laser, and odometry timestamps relative to the IMU, respectively.

[0234] S7.2 Objective Function:

[0235] ;

[0236] in, It is a kernel function for residuals; Q represents the set of keyframe pairs corresponding to the IMU pre-integration constraints; Q represents the set of slices involved in the geometric prior. Represents a set of depth slices; This represents the set of keyframe pairs representing odometer constraints; The weighting coefficients for the wellbore geometric prior residuals are used to adjust the contribution of the geometric prior residuals to the joint optimization objective function. The weighting coefficients for the cross-modal consistency residual term are used to adjust the contribution of the consistency constraint between the laser point cloud and the structured light point cloud to the joint optimization objective function.

[0237] S7.3, Update and retrieve values:

[0238] S7.3-1 Calculation of the number of residual segments and dimensions:

[0239] Inside the sliding window:

[0240] Number of IMU residual segments Dimension of each residual segment ;

[0241] Geometric factors Dimension of each residual segment ;

[0242] Cross-modal factor number Dimension of each residual segment ;

[0243] definition:

[0244] ;

[0245] ;

[0246] ;

[0247] In the above formula, dividing by the quantity and dimension is used to calculate the unit energy of the three terms within the sliding window, and the calculation is performed once after each sliding window iteration;

[0248] S7.3-2, Weight Calculation:

[0249] Calculate the target weights to make the geometry / cross-modal and IMU of the same order of magnitude:

[0250] ;in, Pick Prevent division by zero;

[0251] S7.4, Parameter Update Output:

[0252] The joint optimization objective function is a nonlinear weighted least squares problem with a robust kernel function. It is solved using a combination of iterative reweighted least squares and the Gauss-Newton method: in each iteration, the residual terms are linearized to the first order, and the equivalent weights for the current residuals are calculated based on the robust kernel function. Each residual term is then transformed into a weighted quadratic form, and a normal equation is constructed to solve for the state increment, thereby reducing the impact of outliers on pose, extrinsic parameters, and time offset estimation.

[0253] S7.4-1, Window State Parameters:

[0254] ;

[0255] S7.4-2, External Parameters:

[0256] .

[0257] Furthermore, the specific process of S8 is as follows:

[0258] S8.1, Local element division:

[0259] After preprocessing the acquired structured light and laser point clouds, the point clouds are divided into two-dimensional local units for local feature calculation:

[0260] S8.1-1, Slicing along the wellbore axis:

[0261] ;

[0262] in, This is the initial value along the axial direction. Number the slices. The axial distance in each slice and The minimum and maximum values ​​along the axial direction for each slice;

[0263] point Projected onto the axis direction:

[0264] ;

[0265] like Then point Belongs to the Layer axial slice; if point If the boundary conditions of any axial slice are not met, the point is determined to be outside the current detection axial range or is considered an invalid point and will not participate in the subsequent local unit feature calculation. If the point is located at the common boundary of two slices, it is assigned to the latter slice according to the principle of left closed and right open to avoid the same point being counted repeatedly in multiple slices.

[0266] S8.1-2, Dividing by circumferential angle:

[0267] Partitioning the m-th layer according to its circumference, first calculate the points. The projection vector in the cross-sectional plane perpendicular to the wellbore axis:

[0268] ;

[0269] ;

[0270] in, Point The projection vector in the plane of the wellbore cross-section after removing the axial component; This represents the unit vector along the wellbore axis. Indicates a reference point on the shaft axis;

[0271] Establish a two-dimensional orthogonal basis in a cross-sectional plane perpendicular to the wellbore axis. and And calculate the projection vector. Coordinates in this two-dimensional basis:

[0272] ;

[0273] ;

[0274] in, and These are mutually orthogonal unit basis vectors within the cross-sectional plane; and Representing points respectively Two-dimensional coordinate components within the plane of the wellbore cross-section;

[0275] point The circumferential angle is defined as:

[0276] ;

[0277] in, Point The corresponding circumferential angle within the cross-section of the wellbore. Represent the arctangent function in the four quadrants; if the calculated value is... Then let:

[0278] ;

[0279] Let the number of circumferential divisions be Then the angular width of each circumferential partition is:

[0280] ;

[0281] No. The starting angle of each circumferential partition is defined as:

[0282] ;

[0283] in, Indicates the first The starting angle of each circumferential partition. This indicates the angular width of each circumferential partition. Indicates the total number of circumferential partitions;

[0284] Therefore, the first The first axial slice and the first Each circumferential partition together forms a two-dimensional local unit. For structured light point clouds and laser point clouds, the corresponding local units are defined as follows:

[0285] ;

[0286] ;

[0287] in For structured light point cloud collection, A collection of laser point clouds; and The first The first axial slice, the first The structured light point cloud within a circumferential partition is a local unit of the laser point cloud; and These are the i-th points in the structured light point cloud and the laser point cloud, respectively; Represents the projected coordinates of a point along the axial direction of the wellbore; This indicates the circumferential angle of a point within the cross-section of the well shaft.

[0288] S8.2 Crack detection: Utilizing the high density of structured light point clouds and their sensitivity to fine surface undulations, features such as local normal changes, point cloud density changes, and depth abrupt changes are extracted to identify crack areas. S8.2-1 Calculate the local normal vector:

[0289] For structured light point cloud units For each point Use neighborhood Fit the plane and calculate the normal vector:

[0290] ;

[0291] Calculate the normal gradient of the cell:

[0292] ;

[0293] in Let i be the angle between the normals of point i and its neighboring point j;

[0294] S8.2-2, Local deep mutation:

[0295] Calculate the distance from the structured light spot to the center axis of the unit cell. Then calculate the local depth standard deviation:

[0296] ;

[0297] in The average radial distance of the element. The larger the value, the more obvious the crack or surface unevenness;

[0298] S8.2-3, Calculation of point density index:

[0299] ;

[0300] in, For the first The first axial slice, the first The average radius or average radial distance of each circumferential zone corresponds to a local area of ​​the well wall.

[0301] The following is calculated from the average radial distance from a point within this local unit to the wellbore axis:

[0302] ;

[0303] in, Point The radial distance to the center axis of the wellbore; if the point density is low but the normal phase gradient is high, it is a candidate fracture location;

[0304] S8.2-4 Crack Judgment Scoring Formula

[0305] ;

[0306] in, , and These are the weighting coefficients corresponding to the local normal variation feature, radial discrete feature, and point cloud density feature, respectively. They are used to adjust the contribution ratio of each feature in the crack judgment score. They are all real numbers greater than zero and can be normalized. The weighting coefficients can be set based on experience, determined by calibration experimental results, or obtained through sample training.

[0307] Threshold determination:

[0308] ;

[0309] in, This is the crack detection threshold, used to score the crack detection of local detection units. Perform threshold discrimination; when Greater than the crack detection threshold When this happens, the corresponding local detection unit is identified as a crack area;

[0310] S8.3, Leakage detection:

[0311] By leveraging the geometric stability of laser point clouds and their minimal impact from surface reflection and local pores, local geometric stability and curvature features are extracted. Combined with the abnormal responses of structured light point clouds, such as increased porosity, which are prone to occur under wet conditions, the determination of wet areas or areas of leakage is achieved.

[0312] S8.3-1 Calculation of laser geometric stability:

[0313] For laser data unit Calculate the local radial deviation:

[0314] ;

[0315] in, The number of laser point clouds within a local unit. R represents the distance from each local point cloud to the axis, and R is the fitting radius. The leakage area is typically obtained as follows: Smaller;

[0316] Laser stability weight:

[0317] ;

[0318] S8.3-2, Calculation of structured light porosity:

[0319] For structured light data units :

[0320] ;

[0321] in, The number of structured light point clouds within a local unit; high porosity represents potential leakage areas.

[0322] S8.3-3, Calculation of local curvature of laser beam:

[0323] Laser local unit ;

[0324] Each of the points:

[0325] ;

[0326] Calculate the local centroid:

[0327] ;

[0328] Calculate the covariance matrix:

[0329] ;

[0330] Calculate the eigenvalues ​​of the covariance matrix, assuming the eigenvalues ​​satisfy:

[0331] ;

[0332] in, For the minimum change along the normal direction, This refers to the change along the tangent plane.

[0333] Local curvature formula:

[0334] ;

[0335] S8.3-4, Leakage Assessment:

[0336] ;

[0337] in, For the first The first axial slice, the first The leakage judgment score of each circumferential zone corresponds to a local detection unit; it represents the porosity of the structured light point cloud within the local detection unit, which is used to reflect the degree of structured light point cloud loss caused by water film reflection or wet areas. This represents the geometric stability weight of the laser point cloud within the local detection unit, used to characterize the degree to which the laser point cloud maintains geometric continuity under wet conditions; This indicates the normal variation characteristics of the structured light point cloud within the local detection unit, used to characterize the intensity of local surface orientation changes; This indicates the local curvature characteristics of the laser point cloud within the local detection unit, used to characterize the degree of local geometric curvature of the laser point cloud. and These are used to reduce the leakage score when there is excessive local normal change or excessive local curvature, so as to avoid misjudging obvious geometric defects such as cracks, bulges, and misalignments as leakage areas.

[0338] Let the threshold for determining water leakage be... Used to determine and score leakage in local detection units. Perform threshold discrimination; when Greater than the leakage threshold If the condition is met, the corresponding local detection unit is identified as a leaking area; otherwise, it is identified as a non-leaking area or an area requiring further verification. Threshold It can be determined through experience setting, experimental calibration or sample statistical analysis, and can be adjusted according to the testing conditions.

[0339] This invention achieves high-precision 3D reconstruction and defect detection of vertical shaft walls by integrating multi-source information from structured light, laser, and IMU, through time alignment, construction of well wall geometric prior factors, construction of cross-modal consistency factors, construction of IMU pre-integration factors, construction of slowly varying extrinsic parameter prior factors, construction of odometer factors, sliding window joint optimization, and parallel detection of multi-modal defects from point cloud data. Even under complex conditions such as vibration, water film reflection, and dust obstruction, this invention maintains stable spatial relationships between multiple sensors and consistent point cloud stitching, achieving stable operation while balancing accuracy and robustness. It significantly suppresses extrinsic parameter drift and point cloud distortion, solving the measurement discontinuity and deformation assessment distortion problems caused by data quality fluctuations in traditional single or relative stitching schemes, providing technical support for accurate acquisition of well wall morphology and defect detection. Attached Figure Description

[0340] Figure 1 This is a flowchart of the method of the present invention;

[0341] Figure 2 This is a flowchart of the process of S8 of the present invention. Detailed Implementation

[0342] The invention will now be further described with reference to the accompanying drawings.

[0343] like Figure 1 As shown, a method for detecting the health of a vertical shaft wall based on multi-source data fusion includes the following steps:

[0344] S1. Time alignment: Establish a unified time reference, map the timestamps of the laser scanner, structured light camera and odometer to the unified time reference, and interpolate the mapped observations to obtain synchronous observations on the same time axis as the IMU status.

[0345] S2. Constructing prior geometric factors for the wellbore: Based on the spatial distribution of laser point clouds in the wellbore reference coordinate system, fit the geometric model of the cylindrical surface of the wellbore and construct the geometric prior residuals from the points to the geometric model;

[0346] S3. Construct cross-modal consistency factor: After transforming the laser point cloud and structured light point cloud to a unified coordinate system, slice or layer them according to the axial depth of the well shaft, and construct cross-modal consistency residuals for the cross-sectional parameters of the two modal point clouds within the same depth slice.

[0347] S4. Constructing IMU pre-integration factors: Based on the angular velocity and acceleration measurements of the IMU between adjacent keyframes, pre-integration is performed to construct the kinematic residuals between the states of adjacent keyframes.

[0348] S5. Constructing Slowly Varying Extrinsic Parameter Prior Factors: Set slowly varying prior residuals for the extrinsic parameters between the laser scanner and the IMU, and between the structured light camera and the IMU, to constrain the variation of the extrinsic parameters within the sliding window;

[0349] S6. Constructing the odometer factor: Based on the absolute rope length output by the odometer, construct the axial displacement residual to provide a dimensional constraint on the axial depth of the wellbore.

[0350] S7. Sliding window joint optimization: The pose, velocity, IMU zero bias, extrinsic parameters, and time offset within the window are the variables to be estimated. The residuals from steps S2 to S6 and the marginalized prior residuals are combined to form the objective function and iteratively minimized to obtain the updated pose, extrinsic parameters, and time offset.

[0351] S8, such as Figure 2 As shown, the method for parallel detection of multimodal defects using point cloud data utilizes the high density and sensitivity of structured light point clouds to fine surface undulations to detect crack defects; it also utilizes the geometric stability of laser point clouds and their minimal impact from surface reflection and local pores, combined with the abnormal responses of structured light point clouds, such as increased porosity, which are prone to occur under wet conditions, to determine wet or leaking areas.

Claims

1. A method for detecting the health of a vertical shaft wall based on multi-source data fusion, characterized in that, The steps include the following: S1. Time alignment: Establish a unified time reference, map the timestamps of the laser scanner, structured light camera and odometer to the unified time reference, and interpolate the mapped observations to obtain synchronous observations on the same time axis as the IMU status. S2. Constructing prior geometric factors for the wellbore: Based on the spatial distribution of laser point clouds in the wellbore reference coordinate system, fit the geometric model of the cylindrical surface of the wellbore and construct the geometric prior residuals from the points to the geometric model; S3. Construct cross-modal consistency factor: After transforming the laser point cloud and structured light point cloud to a unified coordinate system, slice or layer them according to the axial depth of the well shaft, and construct cross-modal consistency residuals for the cross-sectional parameters of the two modal point clouds within the same depth slice. S4. Constructing IMU pre-integration factors: Based on the angular velocity and acceleration measurements of the IMU between adjacent keyframes, pre-integration is performed to construct the kinematic residuals between the states of adjacent keyframes. S5. Constructing Slowly Varying Prior Factors for Extrinsic Parameters: Set slowly varying prior residuals for the extrinsic parameters between the laser scanner and the IMU, and between the structured light camera and the IMU, to constrain the variation of the extrinsic parameters within the sliding window; S6. Constructing the odometer factor: Based on the absolute rope length output by the odometer, construct the axial displacement residual to provide a dimensional constraint on the axial depth of the wellbore. S7. Sliding window joint optimization: The pose, velocity, IMU zero bias, extrinsic parameters, and time offset within the window are the variables to be estimated. The residuals from steps S2 to S6 and the marginalized prior residuals are combined to form the objective function and iteratively minimized to obtain the updated pose, extrinsic parameters, and time offset. S8. Parallel detection of multimodal defects using point cloud data: Utilizing the high density and sensitivity to fine surface undulations of structured light point clouds for crack detection; leveraging the geometric stability of laser point clouds and their minimal impact from surface reflection and localized pores, combined with the abnormal responses of structured light point clouds, such as increased porosity under humid conditions, to determine humid or leaking areas.

2. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S1 is as follows: S1.1 Unified Clock: A unified clock source is used for timestamp calibration of the laser scanner, structured light camera, inertial measurement unit, and odometer. The unified clock source is the same master control clock distribution + hardware trigger input, ensuring that the data of the laser scanner, structured light camera, IMU and odometer have the same time reference. S1.2 Software Time Alignment: Based on a unified clock, using the IMU as a high-frequency reference, the observation timestamps of the structured light camera, laser scanner, and odometry are aligned, mapping the observations to the time intervals between two adjacent IMU frames. The IMU states are then interpolated to obtain the pose at the corresponding time for subsequent residual construction. The process is described as follows: ； ； ； in, These are the timestamps from the structured light camera, laser scanner, and odometer, respectively. These represent the time deviations of the observation times of the structured light camera, laser scanner, and odometer relative to the IMU clock. This indicates the corresponding time point on the IMU timeline for observations by the structured light camera. This indicates the corresponding moment on the IMU timeline observed by the laser scanner. This indicates the corresponding moment on the IMU time axis for the odometer observation.

3. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S2 is as follows: S2.1 Constructing the cylindrical surface model: S2.1-1, in the window Within, for each keyframe Take the corresponding laser point Then, use the pose and extrinsic parameters of that keyframe to transform the laser point into the world frame: ； in, Let the laser point be located in world coordinates at the k-th keyframe within the window. for The pose transformation from the body coordinate system to the world coordinate system at any given moment. This is an extrinsic parameter transformation from the laser coordinate system to the body coordinate system. Let be the laser point in the laser coordinate system at the k-th keyframe; The resulting window point set is obtained by aggregation: ； in, This refers to all laser points within window W in the world coordinate system. S2.1-2. Layer according to the z-axis to obtain a set of slices. : ； in, Number the slices. Let z be the z-coordinate of the laser point in the world coordinate system within window W. The minimum z-coordinate value of the laser point in the world coordinate system within window W. The thickness of the slice. This is the m-th slice; S2.1-3. In each slice, project the point cloud onto a plane perpendicular to the axis and fit a circle using the RANSAC algorithm: output the center point of each slice. , The z-coordinate is the average value of the z-coordinate in this slice, and the radius of the circle is... Perform axis fitting and set the center point of each layer. Perform linear fitting: ； The axis to be fitted is represented as: ； in, For a point on the fitted axis, It is a direction vector. The scale value; point The square of the shortest distance to the line is: ； Perform least-squares line fitting: ； The average of the fitted radii for each layer is used as the radius of the window. Finally, the cylindrical parameters of the window are obtained. ; S2.2 Introduction of slice quality gating: Given that downhole dust, water film reflection, and occlusion factors cause sparse laser point clouds or the presence of anomalous echoes, to avoid unreliable fitting misleading joint optimization, slice quality gating is introduced for cylindrical fitting, forming a slice set participating in geometric priors. ; Calculate the fit quality index for each slice: S2.2-1, Interior Point Ratio: ； in, This represents the number of interior points when fitting a circle using the RANSAC algorithm. This represents the number of point clouds within the slice; , Take a value of 0.5 to 0.8; S2.2-2, Coverage: Project the slice points onto the cross-sectional plane, calculate the polar angle distribution of the points, and obtain the coverage angle. : ； in, , Take a value of 0.3 to 0.5; S2.2-3, Fitting Residual: The root mean square error of the fit for all points is: ； ； in, For this slice of point cloud set, For a point in this slice set, Let be a point on the axis. This represents the radial distance residual from the point to the cylindrical surface; , For laser ranging noise, Values ​​range from 2 to 5; S2.2-4, when satisfied When, the slice is included If the geometric prior factors are not constructed, then the geometric prior factors for that slice are constructed. S2.3 Computational geometric prior residuals: cylinder parameters For any point The axis is formed by Definition, point Shortest distance to the axis: ； ； in, , The parameters of the cylinder model obtained by fitting within the window include the cylinder axis. and cylinder radius , Point Relative to cylindrical model parameters The geometric prior residuals; S2.4, For each Take the point cloud set within the slice. Constructing geometric prior residuals: ； in, For the corresponding Point cloud within a slice.

4. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S3 is as follows: S3.1 Calculate the surface-to-surface residual vector: within the depth range covered by the corresponding window. First, the laser point cloud and structured light point cloud are preprocessed by denoising, filtering, and coordinate unification. Then, the point cloud is divided into several depth slices according to the shaft axial depth, forming a depth slice set: ； in, Let j represent the j-th depth slice, and n represent the number of depth slices; For slice unit Point set within The parameters are as follows: S3.1-1 Calculate the centroid: ； Where N is the number of point clouds involved in the calculation in the corresponding slice; S3.1-2, Calculate the covariance: ； S3.1-3, Perform eigenvalue decomposition : Eigenvector corresponding to the smallest eigenvalue That is, the plane normal: ； in, Let be the i-th eigenvalue of the covariance matrix S. For eigenvalues The corresponding feature vector; S3.1-4, Planar Offset: ； Planar offset difference directly reflects the distance difference along the normal direction. Used to constrain translational differences along the normal direction; S3.1-5, Normal Consistency Residual : ； in, Let be the unit normal vector of the locally fitted plane of the laser point cloud. Let be the unit normal vector of the locally fitted plane of the structured light point cloud. This represents the cross product operation of vectors; if the two normals are parallel, the cross product result is 0, i.e. It is a three-dimensional vector residual; S3.1-6, Consistent planar position : The same plane can be written as or In comparison First, ensure that the two normal directions are consistent; otherwise... The inversion of the sign will cause the residual to jump. and Let be the plane offsets of the local fitting planes of the laser point cloud and the structured light point cloud, respectively, and let: ； If not satisfied, then: ； ； S3.1-7, Surface-to-Surface Residual Vector: ； in, This represents the cross-modal surface-to-surface consistency residual vector between the laser point cloud and the structured light point cloud within the corresponding slice z; S3.2, Weight the residuals with confidence level: The problem of residual vectors automatically weakening when data differences exist is addressed by calculating confidence weights based on the cross-sectional quality of structured light point clouds and laser point clouds within each depth slice. This automatically reduces the contribution of slices with poor data quality during the optimization process. ； in, Let the cross-membrane state confidence weight be the value at depth j. The cross-sectional confidence level of the structured light point cloud within the j-th slice; Let the cross-sectional confidence level of the laser point cloud within the j-th slice be denoted as . Based on normal consistency residuals Consistency residuals with planar position The combination yields the first Cross-modal consistency residual vectors for depth slices: ； in, Indicates the first The cross-modal consistency residual vector between the laser point cloud and the structured light point cloud within a depth slice; This represents the normal consistency residual of the local fitting plane of the two modal point clouds within the slice; This indicates the planar position consistency residual of the local fitting plane of the two modal point clouds within the slice; Obtaining transmembrane consistent residuals Then, it is weighted and represented as: ； in, The transmembrane consistency weighted residual corresponding to the j-th depth slice; Indicates the first The cross-modal consistency residual vector between the laser point cloud and the structured light point cloud within a depth slice; For robust functions, used to reduce the impact of abnormal residuals on the results; The credibility formula is: ； For each slice Projecting points onto the cross-sectional plane, fitting circle parameters to the structured light slices. Fitting circle parameters to laser slices The fitting method uses RANSAC + least squares, where: (1) The in-point ratio when fitting the cross section in RANSAC: ; Let be the number of interior points. This represents the total number of point clouds; use Functions that avoid the influence of extreme values ​​and are mapped to [0,1]: ； in, This indicates the quality index to be normalized. This indicates the lower confidence limit of the quality indicator, when... When the normalization result is zero, the result is taken as 0. This indicates the upper limit of confidence for the quality indicator, when When the normalization result is 1, the result is taken as 1. This means that the calculation result will be limited to the range [0,1]. (2) The root mean square distance from the point to the fitted circle after cross-section fitting is: ; It is the radial error from the point to the cross section; (3) Parameters This is a scaling parameter, taken as a multiple of the sensor's ranging standard deviation, to reflect the rapid decay of reliability when the fitting residual exceeds the noise level. Laser point cloud: ; Institutional Light: ; in, This represents the standard deviation of the ranging noise of the laser scanner. The standard deviation of the ranging noise of the structured light camera; (4) The degree of coverage indicates the extent of coverage of the slice on the cross section, i.e., the angular coverage. This is achieved by projecting the slice points onto the cross section plane and calculating the polar angular distribution of the points to obtain the coverage angle. : ； use To avoid the influence of extreme values, map to [0,1]: ； in, This represents the angle coverage to be normalized. This represents the lower bound of the coverage confidence level, when At that time, it was considered that the slice coverage was insufficient, and the normalization result was set to 0; Indicates the upper limit of coverage confidence, when At that time, the slice was considered to have sufficient coverage, and the normalization result was set to 1; This means that the value is restricted to the range [0,1].

5. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S4 is as follows: S4.1 Determine the parameters to be estimated; S4.1-1, Window Status: ； in, These represent IMU's position in the world system. The following posture, position, and speed; These represent zero bias of the IMU gyroscope and zero bias of the accelerometer, respectively. S4.1-2, Global parameters to be estimated: ； in, These represent the extrinsic parameters from the IMU to the camera and the extrinsic parameters from the IMU to the laser, respectively. These represent the offsets of the camera, laser, and odometry timestamps relative to the IMU, respectively. Connecting adjacent states within the window provides constraints, while zero bias is introduced into the pre-integral model for estimation. ; The IMU pre-integration residuals include rotational residuals, velocity residuals, position residuals, gyro zero-bias residuals, and accelerometer zero-bias residuals. ； S4.2 Parameter Definition: Select a set of keyframes within the sliding window. For each adjacent keyframe An IMU pre-integration factor is constructed for each frame to constrain the relative motion between adjacent keyframes: S4.2-1, Set adjacent keyframes The corresponding timestamps are respectively Since the IMU sampling frequency is usually higher than the keyframe frequency, there are multiple IMU sampling times between two adjacent keyframes, denoted as: ； in, Keyframe IMU sampling times between; This is expressed as the number of IMU samples within this adjacent keyframe interval, and satisfies: ； If the keyframe timestamp does not coincide with the IMU sampling time, time interpolation is used to obtain the IMU status or measurement value at the endpoint; Initialize the pre-integral quantity. ； in, For initialization from keyframes arrive Pre-integral relative rotation; Represents a third-order identity matrix; This is the pre-integral relative velocity increment during initialization; This represents the relative position increment during initialization; 0 indicates the zero vector. S4.2-2, Keyframes arrive At the m-th IMU sampling time, the gyroscope angular velocity measurement and accelerometer measurement values ​​output by the IMU are denoted as follows: and Remove the zero bias from it: ； in, Indicates the first The angular velocity vector after removing the gyroscope zero bias at the m-th IMU sampling time within a keyframe interval; This represents the acceleration vector after removing the accelerometer zero bias at the m-th IMU sampling time within the i-th keyframe interval; This indicates that the gyroscope has zero bias. This indicates that the accelerometer has zero bias. The vector representation of the angular velocity after zero bias is: ； in, , , These represent the deflection angular velocity at the m-th IMU sampling time within the i-th keyframe interval. , , Components in three directions; S4.2-3, Rotational pre-integral update: For each IMU sampling interval Perform rotational pre-integral update: ； The time interval between two adjacent IMU sampling times is defined as: ； in, Indicates from keyframe The cumulative pre-integral relative rotation up to the m-th IMU sampling time; express The exponential mapping on the top is used to convert rotation vectors into rotation matrices; during initialization... As the IMU sampling data is gradually integrated, This indicates that the cumulative relative rotation is no longer equivalent to the identity matrix; S4.2-4, Velocity pre-integral update is as follows: ； in, To extract from keyframes The cumulative relative velocity increment obtained up to the m-th IMU sampling time; at this time Used to transfer a vector from the current time-space system Switch to the initial frame machine system ; S4.2-5, Position pre-integration updated as follows: ； in, Indicates from keyframe The pre-integrated relative position increment accumulated up to the m-th IMU sampling time; This represents the debiased acceleration under the sampling time of the m-th IMU. Transition to the starting keyframe Under the machine system, thus accumulating velocity and position pre-integral quantities in the same coordinate system; Complete keyframes and Integrating over all IMU sampling intervals, we obtain the IMU pre-integration between adjacent keyframes: ； ； ； in, , and These represent keyframes. To keyframe The relative rotation, relative velocity increments, and relative position increments are obtained by pre-integration of the IMU sampled data. S4.3 Residual calculation in: S4.3-1 Rotational residuals in: ； ， ； It is a logarithmic mapping of the rotation matrix, outputting a 3D rotation vector; This represents the relative rotation prediction obtained by integrating the gyroscope measurements after zero biasing from i to i+1. This represents the actual rotation obtained by integrating the gyroscope measurements after removing the zero bias during the period from i to i+1; S4.3-2 Velocity residual in the middle: ； in, and These are represented as keyframes. and The velocity of the IMU in the world coordinate system at any given moment; Represented as the gravitational acceleration vector in the world coordinate system; Represented as keyframe and The time interval between; This indicates that the keyframe is estimated from the current state and then transformed. The relative velocity increment under the machine system; This represents the relative velocity increment obtained from IMU pre-integration; S4.3-3 Positional residuals in: ； in, and These are represented as keyframes. and The position of the IMU in the world coordinate system at any given moment; Keyframe The velocity of the IMU in the world coordinate system at any given moment; This is represented as an estimate obtained from the current state and transformed into a keyframe. The relative position increment under the machine system; Represented as the relative position increment obtained from IMU pre-integration; S4.3-4 Zero-biased residuals in: ， ； in, The gyroscope has zero bias residual; The accelerometer has zero bias residual; and These are represented as keyframes. and The gyroscope is at zero bias at any given moment; and These are represented as keyframes. and The accelerometer bias is zero at any given moment; S4.4, IMU residual vector: Construct adjacent keyframe pairs based on keyframes within the sliding window. : ； in, N-1 is the set of adjacent keyframe pairs that participate in the IMU pre-integration constraint within the sliding window; N-1 is the number of keyframes within the sliding window. For any pair of adjacent keyframes ,in By combining the rotational residual, velocity residual, position residual, gyroscope zero-bias residual, and accelerometer zero-bias residual, the IMU pre-integration residual vector is obtained: 。 6. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S5 is as follows: S5.1, Marginalized Prior Residuals: Marginalized Prior Residuals As the sliding window moves forward, the oldest state variable that moves out of the window is marginalized. Its constraint on the objective function is compressed into a priori factors and applied to the retained variables within the window, forming marginalized prior residuals, expressed as: ； in, To preserve the increment of variables within the window, the matrix with vector Obtained by marginalization elimination, the cost functions before and after marginalization are equivalent in the retained variable space; S5.2, Slowly Varying Extrinsic Parameter Prior Residuals: Maintaining the Keyframe Set within a Sliding Window The body to laser external parameters Define the set of adjacent keyframe pairs within the window: ； For any adjacent keyframe pair Define the extrinsic parameter as a slowly varying prior residual: ； in, The rotation matrix represents the attitude of the laser coordinate system relative to the body coordinate system. The translation vector represents the position of the origin of the laser coordinate system within the body coordinate system. Exterior parameters from the body to the structured light camera ,definition: ； in, for Logarithmic mapping is used to represent rotational differences as three-dimensional vectors; The extrinsic parameter estimates at time i are combined to obtain the saturated prior of the extrinsic parameters: ； in, These are the covariance matrices for the changes in laser extrinsic parameters and the covariance matrices for the changes in structured light extrinsic parameters, respectively, used to allow for slow drift of extrinsic parameters and suppress non-physical jumps; S5.3, Prior Residual Term It is defined as a combination of marginalized prior residuals and slowly varying extrinsic prior residuals. ； in, To balance the relative strengths of the two types of priors, the extrinsic slowly varying prior weight coefficients are used. Reliability of cross-modal constraints Related settings; credibility Credibility of cross-modal sections of slices at various depths within the sliding window In summary, the slowly varying prior weights of the extrinsic parameters are adaptively adjusted according to the following relationship: ； in, As the benchmark weight, This is an adjustment coefficient used to control the magnitude of the increase in prior weights when the confidence level decreases.

7. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S6 is as follows: S6.1 Calculate the odometer increment: The odometer outputs the absolute distance for any IMU at critical moments. The odometer increment is defined as: ； in, This represents the axial displacement increment obtained from the odometer between keyframe i and keyframe j. This represents the absolute distance reading of the odometer at time t; It represents the time offset of the odometer timestamp relative to the IMU timestamp, and is used to map critical moments of the IMU to the odometer time axis; When required When the odometer sampling time does not coincide with the odometer sampling time, linear interpolation is used to obtain the odometer reading at that time. ,but: ； in, and The sampling times of two adjacent odometers; Indicates the odometer at time The absolute distance reading, Indicates the odometer at time The absolute distance reading, Indicates the time obtained through linear interpolation. The corresponding absolute distance reading from the odometer; S6.2 Calculate the axial displacement increment estimated by the IMU. : Based on the IMU pose obtained through sliding window optimization, calculate the displacement increment along the well shaft axis between two IMU keyframe moments: set up and They represent and The position vector of the IMU system origin in the world coordinate system at any given time. Let represent a predefined unit vector along the wellbore axis. Then, the axial displacement increment estimated by the IMU is: ； in, This represents the unique increment of the wellbore axis obtained by the pose estimation of the IMU between keyframe i and keyframe j; They are respectively Displacement estimated by the IMU at time step; The defined unit vector along the wellbore axis; This indicates that the three-dimensional displacement between two keyframes is projected onto the shaft axial direction. S6.3 Calculate the odometer residual: Define the set of adjacent keyframe pairs within the window: ； For any adjacent keyframe pair Define the extrinsic parameter as a slowly varying prior residual: 。 8. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S7 is as follows: S7.1, Parameters to be estimated: Within the sliding window, IMU pose, velocity, position, IMU bias, sensor extrinsic parameters, and the time offset of each sensor relative to the IMU are used as joint parameters to be estimated. Through IMU motion constraints, wellbore geometry priors, cross-modal consistency constraints, slowly varying extrinsic parameter priors, and odometry residual constraints, IMU pose optimization and extrinsic parameter updates are performed within the same optimization framework. The parameters to be estimated consist of two parts: S7.1-1, Window Status: ； in, These represent IMU's position in the world system. The following posture, position, and speed; These represent zero bias of the IMU gyroscope and zero bias of the accelerometer, respectively. S7.1-2, Global parameters to be estimated: ； in, These represent the extrinsic parameters from the IMU to the camera and the extrinsic parameters from the IMU to the laser, respectively. These represent the offsets of the camera, laser, and odometry timestamps relative to the IMU, respectively. S7.2 Objective Function: ； in, It is a kernel function for residuals; Q represents the set of keyframe pairs corresponding to the IMU pre-integration constraints; Q represents the set of slices involved in the geometric prior. Represents a set of depth slices; The set of keyframe pairs representing odometer constraints; The weighting coefficients for the wellbore geometric prior residuals are used to adjust the contribution of the geometric prior residuals to the joint optimization objective function. The weighting coefficients for the cross-modal consistency residual term are used to adjust the contribution of the consistency constraint between the laser point cloud and the structured light point cloud to the joint optimization objective function. S7.3, Update and retrieve values: S7.3-1 Calculation of the number of residual segments and dimensions: Inside the sliding window: Number of IMU residual segments Dimension of each residual segment ; Geometric factors Dimension of each residual segment ; Cross-modal factor number Dimension of each residual segment ; definition: ； ； ； In the above formula, dividing by the quantity and dimension is used to calculate the unit energy of the three terms within the sliding window, and the calculation is performed once after each sliding window iteration; S7.3-2, Weight Calculation: Calculate the target weights to make the geometry / cross-modal and IMU of the same order of magnitude: ;in, Pick Prevent division by zero; S7.4, Parameter Update Output: The joint optimization objective function is a nonlinear weighted least squares problem with a robust kernel function. It is solved using a combination of iterative reweighted least squares and the Gauss-Newton method: in each iteration, the residual terms are linearized to the first order, and the equivalent weights for the current residuals are calculated based on the robust kernel function. Each residual term is then transformed into a weighted quadratic form, and a normal equation is constructed to solve for the state increment, thereby reducing the impact of outliers on pose, extrinsic parameters, and time offset estimation. S7.4-1, Window State Parameters: ； S7.4-2, External Parameters: 。 9. The method for detecting the health of a vertical shaft wall based on multi-source data fusion according to claim 1, characterized in that, The specific process of S8 is as follows: S8.1, Local element division: After preprocessing the acquired structured light and laser point clouds, the point clouds are divided into two-dimensional local units for local feature calculation: S8.1-1, Slicing along the wellbore axis: ； in, This is the initial value along the axial direction. Number the slices. The distance along the axial direction in each slice and The minimum and maximum values ​​along the axial direction for each slice; point Projected onto the axis direction: ； like Then point Belongs to the Layer axial slice; if point If the boundary conditions of any axial slice are not met, the point is determined to be outside the current detection axial range or is considered an invalid point and will not participate in the subsequent local unit feature calculation. If the point is located at the common boundary of two slices, it is assigned to the latter slice according to the principle of left closed and right open to avoid the same point being counted repeatedly in multiple slices. S8.1-2, Dividing by circumferential angle: Partitioning the m-th layer according to its circumference, first calculate the points. The projection vector in the cross-sectional plane perpendicular to the wellbore axis: ； ； in, Point The projection vector in the plane of the wellbore cross-section after removing the axial component; This represents the unit vector along the wellbore axis. Indicates a reference point on the shaft axis; Establish a two-dimensional orthogonal basis in a cross-sectional plane perpendicular to the wellbore axis. and And calculate the projection vector. Coordinates in this two-dimensional basis: ； ； in, and These are mutually orthogonal unit basis vectors within the cross-sectional plane; and Representing points respectively Two-dimensional coordinate components within the plane of the wellbore cross-section; point The circumferential angle is defined as: ； in, Point The corresponding circumferential angle within the cross-section of the wellbore. Represent the arctangent function in the four quadrants; if the calculated value is... Then let: ； Let the number of circumferential divisions be Then the angular width of each circumferential partition is: ； No. The starting angle of each circumferential partition is defined as: ； in, Indicates the first The starting angle of each circumferential partition. This indicates the angular width of each circumferential partition. Indicates the total number of circumferential partitions; Therefore, the first The first axial slice and the first Each circumferential partition together forms a two-dimensional local unit. For structured light point clouds and laser point clouds, the corresponding local units are defined as follows: ； ； in, For structured light point cloud collection, A collection of laser point clouds; and The first The first axial slice, the first The structured light point cloud within a circumferential partition is a local unit of the laser point cloud; and These are the i-th points in the structured light point cloud and the laser point cloud, respectively; Represents the projected coordinates of a point along the axial direction of the wellbore; This indicates the circumferential angle of a point within the cross-section of the well shaft. S8.2 Crack detection: Utilizing the high density of structured light point clouds and their sensitivity to fine surface undulations, features such as local normal changes, point cloud density changes, and depth abrupt changes are extracted to identify crack areas. S8.2-1 Calculate the local normal vector: For structured light point cloud units For each point Use neighborhood Fit the plane and calculate the normal vector: ； Calculate the normal gradient of the cell: ； in Let i be the angle between the normals of point i and its neighboring point j; S8.2-2, Local deep mutation: Calculate the distance from the structured light spot to the center axis of the unit cell. Then calculate the local depth standard deviation: ； in The average radial distance of the element. The larger the value, the more obvious the crack or surface unevenness; S8.2-3, Calculation of point density index: ； in, For the first The first axial slice, the first The average radius or average radial distance of each circumferential zone corresponds to a local area of ​​the well wall. The following is calculated from the average radial distance from a point within this local unit to the wellbore axis: ； in, Point The radial distance to the center axis of the wellbore; if the point density is low but the normal phase gradient is high, it is a candidate fracture location; S8.2-4 Crack Judgment Scoring Formula ； in, , and These are the weighting coefficients corresponding to the local normal variation feature, radial discrete feature, and point cloud density feature, respectively, used to adjust the contribution ratio of each feature in the crack judgment score. They are all real numbers greater than zero and can be normalized. Threshold determination: ； in, This is the crack detection threshold, used to score the crack detection of local detection units. Perform threshold discrimination; when Greater than the crack detection threshold When this happens, the corresponding local detection unit is identified as a crack area; S8.3, Leakage detection: S8.3-1 Calculation of laser geometric stability: For laser data unit Calculate the local radial deviation: ; in, The number of laser point clouds within a local unit. R represents the distance from each laser local point cloud to the axis, and R is the fitting radius. This is typically used to obtain the distance from the laser local point cloud to the axis. Smaller; Laser stability weight: ； S8.3-2, Calculation of structured light porosity: For structured light data units : ； in, The number of structured light point clouds within a local unit; high porosity represents potential leakage areas. S8.3-3, Calculation of local curvature of laser beam: Laser local unit ; Each of the points: ； Calculate the local centroid: ； Calculate the covariance matrix: ； Calculate the eigenvalues ​​of the covariance matrix, assuming the eigenvalues ​​satisfy: ； in, For the minimum change along the normal direction, This refers to the change along the tangent plane. Local curvature formula: ； S8.3-4, Leakage Assessment: ； in, For the first The first axial slice, the first The leakage judgment score of each circumferential zone corresponds to a local detection unit; it represents the porosity of the structured light point cloud within the local detection unit, which is used to reflect the degree of structured light point cloud loss caused by water film reflection or wet areas. This represents the geometric stability weight of the laser point cloud within the local detection unit, used to characterize the degree to which the laser point cloud maintains geometric continuity under wet conditions; This indicates the normal variation characteristics of the structured light point cloud within the local detection unit, used to characterize the intensity of local surface orientation changes; This indicates the local curvature characteristics of the laser point cloud within the local detection unit, used to characterize the degree of local geometric curvature of the laser point cloud; and These are used to reduce the leakage score when there is excessive local normal change or excessive local curvature, so as to avoid misjudging obvious geometric defects such as cracks, bulges, and misalignments as leakage areas. Let the threshold for determining water leakage be... Used to determine and score leakage in local detection units. Perform threshold discrimination; when Greater than the leakage threshold If the condition is met, the corresponding local detection unit will be identified as a leaking area; otherwise, it will be identified as a non-leaking area or an area requiring further verification.

Images