Combined inertial navigation and vision collaborative underground pipeline internal anomaly positioning method

Abstract

The invention discloses a combined inertial navigation and vision collaborative underground pipeline internal anomaly positioning method. An active crawling type pipeline robot is used for carrying a camera, a light supplementing lamp, an inertial measurement unit IMU, a wheel type odometer and a geomagnetic sensor to carry out pipeline internal data collection and transmission. Carrying out real-time positioning calculation on the pipeline robot by utilizing a positioning method based on self-adaptive cubic Kalman filtering, and giving coordinate data of the pipeline robot at each moment; an image frame collected by a camera is analyzed by using an unsupervised pipeline interior anomaly detection method based on graph model multi-scale feature fitting, and detection of unknown anomaly in a pipeline is realized. And when an abnormal condition is detected in the pipeline, the current robot positioning coordinate is fed back, and accurate feedback of an abnormal position is realized. Accurate coordinate information is provided for pipeline internal abnormity repair, the damage range is reduced, the maintenance efficiency is improved, the maintenance cost is reduced, and the method has high engineering application value.

Description

Technical Field

[0001] The invention relates to a combined inertial navigation and vision-coordinated method for locating anomalies inside an underground pipeline, and belongs to the field of navigation positioning and computer vision. Background Art

[0002] The construction of underground pipeline projects is a concealed project, and the quality of its acceptance is directly related to the safe operation of underground pipelines. Therefore, how to conduct efficient and intelligent inspections of underground pipelines at multiple stages of laying to ensure the safe operation of underground pipelines is a key issue that needs to be solved in the future.

[0003] Using an underground pipeline robot equipped with a camera to enter the interior of the pipeline for non-destructive visual inspection is currently a mainstream means of realizing the diagnosis of the internal status of underground pipelines. However, due to the poor lighting conditions inside the underground pipeline, and the various and hidden abnormalities such as water seepage, damage, and rust inside the pipeline, it is difficult to detect abnormalities inside the underground pipeline. In addition, since the pipeline is buried deep underground and cannot receive satellite positioning signals, it is difficult for the robot to determine the specific location of the abnormal area after detecting the abnormality using visual information, and it is impossible to provide specific positioning information for the repair of underground pipeline abnormalities. Therefore, it is of great significance and value to study an abnormality detection method that is adaptable to the environment and obtains the positioning coordinates of the abnormal area under the condition of satellite signal denial. Summary of the invention

[0004] Purpose of the invention: The technical problem that the present invention actually aims to solve is to propose a method for locating abnormalities inside underground pipelines by combining inertial navigation with vision, overcoming the influence of the narrow interior of underground pipelines, the absence of satellite signals, and poor lighting conditions, and using pipeline robots to intelligently detect abnormal areas and accurately provide underground positioning information of abnormal areas, thereby providing accurate reference information for pipeline operation and maintenance.

[0005] The above purpose is achieved through the following technical solutions:

[0006] A combined inertial navigation and vision-coordinated method for locating anomalies inside an underground pipeline of the present invention comprises the following steps:

[0007] (1) An active crawling pipeline robot is equipped with a camera, a fill light, an inertial measurement unit (IMU) with a built-in three-axis accelerometer and three-axis gyroscope, a wheel odometer, and a geomagnetic sensor to collect and transmit data inside the pipeline;

[0008] (2) The positioning method based on adaptive cubic Kalman filter is used to perform real-time positioning calculation of the pipeline robot and provide the coordinate data Pos of the pipeline robot at each moment. k ;

[0009] (3) An unsupervised pipeline anomaly detection method based on multi-scale feature fitting of a graphical model is used to analyze the image frame Image(k) collected by the camera to detect abnormal conditions inside the pipeline.

[0010] (4) When an abnormal situation is detected in the pipeline through step (3), the current robot positioning coordinate Pos is fed back k , achieving accurate feedback of abnormal positions.

[0011] Furthermore, step (1) specifically comprises the following steps:

[0012] (11) Before the pipeline robot is put into operation, the sensor needs to be calibrated. First, the accelerometer and gyroscope of the inertial measurement unit (IMU) are bias-corrected, and the initial orientation of the geomagnetic sensor is determined using static measurement.

[0013] (12) After the pipeline robot enters the pipeline, it turns on the camera and fill light, and starts to collect the robot's own positioning data and image data. Assume that the data collected at time k are: The acceleration measurement value of the inertial measurement unit IMU in the body coordinate system is They represent the acceleration measurement values ​​of the pipeline robot in the x-axis, y-axis, and z-axis directions in the body coordinate system at time k, respectively. The superscript T represents the transpose of the vector and the angular velocity measurement value of the inertial measurement unit IMU in the body coordinate system. They represent the angular velocity measurement values ​​of the pipeline robot in the x-axis, y-axis, and z-axis directions in the body coordinate system at time k, and the odometer measurement speed in the body coordinate system δs k Indicates the speed measured by the odometer in the forward or backward direction of the pipeline robot and the geomagnetic measurement value in the body coordinate system Respectively represent the geomagnetic measurement values ​​and gravity accelerometer data of the pipeline robot in the east, north, and up directions at time k The image frame Image(k) acquired by the camera; the subscripts x, y, and z respectively represent the "x, y, and z" axes in the body coordinate system, pointing to the east, north, and up directions respectively; the subscript b represents the body coordinate system of the pipeline robot; it is calculated by measuring the speed of wheel rotation; "k" represents the kth moment, and when k = 0, it represents the initial state value;

[0014] (13) The collected data will be transmitted to the ground workstation for processing via the optical fiber cable at the rear end of the active crawling pipeline robot.

[0015] Furthermore, step (2) specifically comprises the following steps:

[0016] (21) Define the state vector of the pipeline robot at time k as X k , Definition Φ k is the pipeline robot pose, are the pitch, roll and yaw angles of the pipeline robot at time k; δΦ k is the navigation angle error of the robot at time k, in, They represent the pitch, roll and yaw angle errors of the pipeline robot at time k respectively; V k represents the speed of the pipeline robot at time k in the navigation coordinate system, They represent the speed of the pipeline robot in the east, north and upward directions at time k; δV k is the velocity error of the pipeline robot in the navigation coordinate system at time k, They represent the velocity errors of the pipeline robot in the east, north, and up directions at time k, respectively; is the bias error of the angular velocity measurement value of the inertial measurement unit IMU of the pipeline robot in the body coordinate system at time k, They respectively represent the bias errors of the angular velocity measurements of the inertial measurement unit IMU in the x-axis, y-axis, and z-axis directions in the body coordinate system of the pipeline robot at time k; represents the IMU accelerometer bias error of the pipeline robot in the body coordinate system at time k, They represent the IMU accelerometer bias errors in the x-axis, y-axis, and z-axis directions of the pipeline robot in the body coordinate system at time k, respectively;

[0017] Define the state error covariance matrix of the pipeline robot at time k as P k , They respectively represent the variance of the pipeline robot's posture at time k, the velocity variance in the navigation coordinate system, the bias error variance of the angular velocity measurement value of the inertial measurement unit IMU of the pipeline robot in the body coordinate system at time k, and the bias error variance of the IMU accelerometer in the body coordinate system;

[0018] Define the process noise covariance matrix of the pipeline robot at time k as Q k , They respectively represent the variance of the angular velocity measurement value of the inertial measurement unit IMU in the body coordinate system and the variance of the acceleration measurement value of the inertial measurement unit IMU in the body coordinate system;

[0019] Define the measurement noise covariance matrix of the pipeline robot at time k as R k , They respectively represent the variance of the speed measured by the odometer in the body coordinate system, the variance of the geomagnetic measurement value in the body coordinate system, and the variance of the acceleration measurement value of the inertial measurement unit IMU in the body coordinate system;

[0020] Define the coordinates of the pipeline robot at time k as Pos k ,Pos k =[x k ,y k ,z k ] T , x k ,y k 、z k They represent the coordinates of the x, y, and z axes in the body coordinate system, respectively, and the superscript T represents the transpose of the vector;

[0021] definition is the transformation matrix from the body coordinate system to the navigation coordinate system at time k;

[0022] (22) Pipeline robot sensor data initialization: Initialize state vector: Initialize the state error covariance matrix: Initialize the process noise covariance matrix: Initialize the measurement noise covariance matrix: Initialize the pipeline robot positioning coordinates: Pos0 = [x0, y0, z0] T , the transformation matrix of the initial posture from the body coordinate system to the navigation coordinate system is defined as:

[0023] (24) Adaptive cubic Kalman filter is used for state prediction, and the specific calculation steps are as follows:

[0024] From d Draw cube sampling points at the intersection of the dimensional unit sphere and the Cartesian coordinate axes and perform a scale transformation:

[0025]

[0026] Among them, e i is a unit vector, the i-1th element is 1, the rest are 0, ξ i Indicates that from n d The result of scaling after the cube sampling points are extracted from the intersection of the unit sphere and the Cartesian coordinate axis, n d Represents the dimension of the space, 2n d Indicates the total number of cube sampling points extracted;

[0027] Project the cube sampling point into the state space to obtain the cube sampling point X at time k-1 i,k-1∣k-1 :

[0028]

[0029] Among them, P k-1∣k-1 is the covariance matrix of the state estimate at time k-1, is the mean of the state estimate at time k-1;

[0030] Then X i,k-1∣k-1 Substituting into the nonlinear dynamic model, the cube after propagation adopts point That is, the predicted state at time k:

[0031]

[0032] Where f(·) is the discrete-time state prediction equation. The state error propagation equation is discretized by numerical integration to calculate the predicted state of the system from time k-1 to time k. The state error propagation equation of the system is:

[0033]

[0034] in, is the time rate of change of the pitch, roll, and heading attitude errors from time k-1 to time k, is the time rate of change of the velocity error from time k-1 to time k;

[0035] δΦ k-1 × indicates opposition to matrix form and has:

[0036] is the angular velocity of the navigation system at time k-1, is the angular velocity error in the navigation system at time k-1,

[0037] is the angular velocity of the Earth's rotation, is the error of the Earth's rotation angular velocity:

[0038]

[0039] Among them, Ω is the rotation rate of the earth, L is the latitude, and δL is the latitude error;

[0040] is the navigation error angular velocity at time k:

[0041]

[0042] Where R is the mean radius of curvature of the earth; h r is the height of the pipeline robot; are the eastward and northward velocities of the pipeline robot at time k-1 respectively; is the speed error of the pipeline robot in the east and north directions at time k-1; δh r is the height error of the pipeline robot;

[0043] N k-1 is the compensation term caused by the earth's rotation and navigation error at time k-1, M k-1 is the compensation term caused by the inertial system error at time k-1,

[0044] Calculate the predicted state mean:

[0045]

[0046] Compute the forecast error covariance:

[0047]

[0048]

[0049] Process noise covariance matrix Q k-1 Adaptive adjustment based on IMU raw gyroscope data: Q k-1 =Q0×10 ξn , improve the stability of the filter under different conditions, where Q0 is the initial value, is the adaptive algorithm parameter, is the original angular velocity measured by the gyroscope;

[0050] (24) Adaptive cubic Kalman filtering is used for measurement update, and the specific calculation steps are as follows:

[0051] (a) Calculate the pitch angle θ of the pipeline robot at time k by calculating the accelerometer and geomagnetic sensor data k , yaw angle ψ k , roll angle β k :

[0052]

[0053] Where μ is the relative magnetic permeability, They represent the acceleration components in the x, y, and z directions of the body coordinate system measured by the IMU accelerometer at time k, is the projection of the geomagnetic field at time k in the coordinate system of the pipeline robot, |B| is the total intensity of the geomagnetic field;

[0054] Calculate the transformation matrix from the body coordinate system to the navigation coordinate system at time k

[0055]

[0056] (b) Setting the threshold ξ g , determine whether the current IMU data is available:

[0057]

[0058] Among them, g z is the acceleration due to gravity, × represents the cross product, ξ g is a constant threshold;

[0059] (c) Obtain the observation vector based on the odometer, accelerometer, and magnetometer sensors: The speed measurement error is is the speed in the navigation coordinate system at time k, is the transformation matrix from the fuselage coordinate system to the navigation coordinate system at time k, is the speed in the fuselage coordinate system measured by the odometer at time k, is the projection of the geomagnetic field at time k in the coordinate system of the pipeline robot;

[0060] (d) Calculate the new cube sampling point X based on the predicted state covariance i,k∣k-1 :

[0061]

[0062] Among them, X i,k∣k-1 represents the i-th cubic point, represents the cubic point after the propagation of the predicted state covariance at time k-1, P k∣k-1 is the predicted state covariance matrix, which represents the uncertainty of the state variables, E{X k∣k-1} is the predicted state mean, which is used to adjust the center position of the cube point;

[0063] Calculate the cube sampling points P in the measurement space using the measurement model xz,k∣k-1 :

[0064] Z i,k∣k-1 =h(X i,k∣k-1 ) (13)

[0065] Where h(·) is the measurement model;

[0066] Calculates predicted measurement estimates

[0067]

[0068] Calculate the covariance matrix S for the uncertainty between the predicted and actual measurementsk∣k-1 :

[0069]

[0070] Among them, R k is the measurement noise covariance matrix;

[0071] Calculate the state and measurement cross covariance matrix P xz,k∣k-1 :

[0072]

[0073] Calculate the Kalman gain K k :

[0074]

[0075] (25) Correct the pipeline robot status error:

[0076] (a) Update the state estimate and calculate the predicted state mean E{X k|k}:

[0077]

[0078] Among them, Z k is the actual measurement value of the sensor, is the predicted measurement;

[0079] (b) Update the error covariance P k|k :

[0080] P k|k =P k∣k-1 -K k P zz,k∣k-1 (K k ) T (19)

[0081] Among them, P k∣k-1 is the predicted state covariance matrix, P zz,k∣k-1 Represents the forecast measurement error covariance matrix:

[0082]

[0083] (c) Establish the velocity error correction equation and attitude error correction equation, and update the system state parameters:

[0084] Establish the speed error correction equation:

[0085]

[0086] in, is the corrected velocity value at time k, is the speed value calculated at time k, [X k (5),X k (6),X k (7)] is the velocity error correction obtained by Kalman filtering at time k;

[0087] Establish the attitude error correction equation:

[0088]

[0089] in, is the modified k-time attitude transformation matrix, is the original k-1 moment posture transformation matrix, [X k (1),X k (2),X k (3)] ​​is the attitude error correction obtained by Kalman filtering;

[0090] (26) Update the location information of the pipeline robot:

[0091]

[0092] Among them, Pos k-1 represents the position information of the pipeline robot at time k-1, that is, the position at the previous moment, and Δt represents the time difference between time k and time k-1.

[0093] Furthermore, step (3) specifically comprises the following steps:

[0094] (33) Load the pre-trained ResNet-50 network model as a feature extractor;

[0095] (34) Construct the feature representation of the normal image set, that is, the gallery set F G ;

[0096] (a) Load a set of M defect-free underground pipeline interior images as the normal image set D = {D1, D2, D3, ..., D M}, each image Among them, 1≤i≤M represents the total number of normal images, w, h, c are the width, height and length of the image respectively;

[0097] (b) ResNet-50 network is used to extract the feature map of the normal image set. Suppose the feature map extracted by the lth layer of the network is Among them, l∈{1,2,3,4}, w l ,h l Respectively represent the width and height of the feature map, c l Represents the number of channels of the feature map of the lth layer. The feature vector of each pixel point p on the lth layer is expressed as 1≤p≤p l , where p l =w l ×h l Represents the total number of feature points in the feature map of layer l;

[0098] (c) For all normal images without defects D i ,1≤i≤M and all layers l∈{1,2,3,4}, summarize the feature vectors of all feature points to form a local gallery set

[0099]

[0100] (d) Perform global average pooling on each image to obtain the global gallery set

[0101]

[0102] (e) Combine the local gallery set with the global gallery set to construct the final complete gallery set

[0103] (33) Use the ResNet-50 network to collect images in real time from the pipeline robot camera Perform feature extraction, w, h, c are the width, height and length of the image respectively, and construct the query set F Q ;

[0104] (a) Using ResNet-50 network to extract real-time acquired images The feature map of the network is as follows: Among them, l∈{1,2,3,4}, w l ,h l Respectively represent the width and height of the feature map, c l Represents the number of channels of the feature map of the lth layer; the feature vector of each pixel point p on the lth layer is expressed as 1≤p≤p l , where p l =w l ×h l Represents the total number of feature points in the feature map of layer l;

[0105] (b) For Image(k) and all layers l∈{1,2,3,4}, the feature vectors of all feature points are aggregated to form a local query set

[0106]

[0107] (c) Perform global average pooling on each image to obtain the global query set

[0108]

[0109] (d) Combine the local query set with the global query set to construct the final complete query set

[0110] (34) For each query point y i ∈F Q , in Gallery Set F G Find k s Nearest neighbors {f i}, and calculate the cosine similarity a i,j :

[0111]

[0112] Select the first K largest similarity f j As y i The nearest neighbor, that is, N i = {v g,1 ,v g,2 ,…,v g,j ,…,v g,K}, the vertices corresponding to these nearest neighbor features are represented as v g,j ∈V G , V G is a set of graph vertices, and constructs an adjacency matrix A, with a i,j ∈A:

[0113]

[0114] (35) For query point y i , using K nearest neighbors for weighted feature fitting:

[0115]

[0116] Calculate the anomaly score s for each vertex based on the cosine distance i :

[0117]

[0118] Get the anomaly score matrix And use bilinear interpolation method to Zoom in to the same size as the original test image:

[0119]

[0120] According to the weights of different layers, the weighted multi-scale anomaly score matrix Ω(Image(k)) is calculated:

[0121]

[0122] in, Represents the weight coefficient of the feature map of the lth layer;

[0123] (37) Set the threshold τ1 to determine whether the pixel is abnormal. If Ω(Image(k),p)>τ1, the pixel p is marked as an abnormal pixel. Set the threshold τ2. If the maximum abnormal value of the entire image max(Ω(Image(k)))>τ2, Image(k) is judged as an abnormal image, and a binary abnormal segmentation map Defect(Image(k)) is generated:

[0124]

[0125] Compared with the prior art, the present invention has the following beneficial effects:

[0126] (1) The present invention utilizes an active crawling pipeline robot equipped with a visible light camera, an IMU, and a coding wheel to obtain positioning information and visual information inside the pipeline, which can effectively replace manual acquisition of long-distance pipeline internal information and solve the problem of unclear internal pipeline information.

[0127] (2) The present invention proposes a combined inertial navigation positioning method that combines IMU, encoder wheel, geomagnetic signal and gravimeter, introduces geomagnetic and gravity sensors to obtain the observation vector of the filter, and uses an improved spherical radial integral Kalman filter (ACKF) method to solve the nonlinear problem of the filter attitude error. By using a low-cost, miniaturized, and low-power sensor combination, the crawling trajectory of the robot in the pipeline can be determined to obtain accurate three-dimensional positioning information of the underground pipeline.

[0128] (3) The present invention proposes an unsupervised pipeline anomaly detection method based on multi-scale feature fitting of a graphical model, which can effectively detect and locate different types of and unknown anomalies in the pipeline, and solve the problems of complexity, diversity, strong randomness, and small number of samples of anomalies inside the pipeline.

[0129] (4) By combining the positioning information of the robot during crawling with the results of intelligent visual detection, the accurate positioning information of abnormalities inside underground pipelines that cannot be reached by human power can be determined, providing accurate coordinate information for excavation and repair, reducing the scope of damage, improving maintenance efficiency, and reducing maintenance costs. BRIEF DESCRIPTION OF THE DRAWINGS

[0130] Figure 1 This is a block diagram of the method for locating abnormalities inside underground pipelines using combined inertial navigation and vision collaboration.

[0131] Figure 2 Schematic diagram of the positioning method based on adaptive cubic Kalman filter.

[0132] Figure 3 Schematic diagram of the unsupervised in-pipeline anomaly detection method based on multi-scale feature fitting of the graphical model. DETAILED DESCRIPTION

[0133] The method for locating abnormalities inside underground pipelines by combining inertial navigation and vision in the present invention is as follows: Figure 1 The specific operation process is as follows:

[0134] (1) An active crawling pipeline robot equipped with a camera, fill light, inertial measurement unit (IMU) (with built-in three-axis accelerometer and three-axis gyroscope), wheel odometer, and geomagnetic sensor is used to collect and transmit data inside the pipeline.

[0135] (11) Before the robot is running, the sensors need to be calibrated. First, the accelerometer and gyroscope of the IMU are biased and the initial orientation of the geomagnetic sensor is determined using static measurements.

[0136] (12) After the pipeline robot enters the pipeline, it turns on the camera and fill light, and starts to collect the robot's own positioning data and image data. Assume that the data collected at time k are: The acceleration measurement value of the inertial measurement unit IMU in the body coordinate system is They represent the acceleration measurement values ​​of the pipeline robot in the x-axis, y-axis, and z-axis directions in the body coordinate system at time k, respectively. The superscript T represents the transpose of the vector and the angular velocity measurement value of the inertial measurement unit IMU in the body coordinate system. They represent the angular velocity measurement values ​​of the pipeline robot in the x-axis, y-axis, and z-axis directions in the body coordinate system at time k, and the odometer measurement speed in the body coordinate system δs k Indicates the speed measured by the odometer in the forward or backward direction of the pipeline robot and the geomagnetic measurement value in the body coordinate system Respectively represent the geomagnetic measurement values ​​and gravity accelerometer data of the pipeline robot in the east, north, and up directions at time k The image frame Image(k) acquired by the camera; the subscripts x, y, and z respectively represent the "x, y, and z" axes in the body coordinate system, pointing to the east, north, and up directions respectively; the subscript b represents the body coordinate system of the pipeline robot; it is calculated by measuring the speed of wheel rotation; "k" represents the kth moment, and when k=0, it represents the initial state value.

[0137] (13) The collected data will be transmitted to the ground workstation for processing via the optical fiber cable at the rear end of the active crawling pipeline robot.

[0138] (2) Figure 2 As shown in the figure, the positioning method based on adaptive cubic Kalman filter is used to solve the real-time positioning of the pipeline robot, and the coordinate data Pos of the pipeline robot at each moment is given. k . (twenty one)

[0140] (21) Define the state vector of the pipeline robot at time k as X k , Definition Φ k is the pipeline robot pose, are the pitch, roll and yaw angles of the pipeline robot at time k; δΦ k is the navigation angle error of the robot at time k, in, They represent the pitch, roll and yaw angle errors of the pipeline robot at time k respectively; V k represents the speed of the pipeline robot at time k in the navigation coordinate system, They represent the speed of the pipeline robot in the east, north and upward directions at time k; δV k is the velocity error of the pipeline robot in the navigation coordinate system at time k, They represent the velocity errors of the pipeline robot in the east, north, and up directions at time k, respectively; is the bias error of the angular velocity measurement value of the inertial measurement unit IMU of the pipeline robot in the body coordinate system at time k, They respectively represent the bias errors of the angular velocity measurements of the inertial measurement unit IMU in the x-axis, y-axis, and z-axis directions in the body coordinate system of the pipeline robot at time k; represents the IMU accelerometer bias error of the pipeline robot in the body coordinate system at time k, They represent the IMU accelerometer bias errors in the x-axis, y-axis, and z-axis directions of the pipeline robot in the body coordinate system at time k, respectively;

[0141] Define the state error covariance matrix of the pipeline robot at time k as P k , They respectively represent the variance of the pipeline robot's posture at time k, the velocity variance in the navigation coordinate system, the bias error variance of the angular velocity measurement value of the inertial measurement unit IMU of the pipeline robot in the body coordinate system at time k, and the bias error variance of the IMU accelerometer in the body coordinate system;

[0142] Define the process noise covariance matrix of the pipeline robot at time k as Q k , They respectively represent the variance of the angular velocity measurement value of the inertial measurement unit IMU in the body coordinate system and the variance of the acceleration measurement value of the inertial measurement unit IMU in the body coordinate system;

[0143] Define the measurement noise covariance matrix of the pipeline robot at time k as R k , They respectively represent the variance of the speed measured by the odometer in the body coordinate system, the variance of the geomagnetic measurement value in the body coordinate system, and the variance of the acceleration measurement value of the inertial measurement unit IMU in the body coordinate system;

[0144] Define the coordinates of the pipeline robot at time k as Pos k ,Pos k =[x k ,y k ,z k ] T , x k ,y k 、z k They represent the coordinates of the x, y, and z axes in the body coordinate system, respectively, and the superscript T represents the transpose of the vector;

[0145] definition is the transformation matrix from the body coordinate system to the navigation coordinate system at time k.

[0146] (22) Pipeline robot sensor data initialization: Initialize state vector: Initialize the state error covariance matrix: Initialize the process noise covariance matrix: Initialize the measurement noise covariance matrix: Initialize the pipeline robot positioning coordinates: Pos0 = [x0, y0, z0] T , the transformation matrix of the initial posture from the body coordinate system to the navigation coordinate system is defined as:

[0147] (23) Adaptive cubic Kalman filter is used for state prediction, and the specific calculation steps are as follows:

[0148] From d Draw cube sampling points at the intersection of the dimensional unit sphere and the Cartesian coordinate axes and perform a scale transformation:

[0149]

[0150] Among them, e i is a unit vector, the i-1th element is 1, the rest are 0, ξ i Indicates that from n dThe result of scaling after the cube sampling points are extracted from the intersection of the unit sphere and the Cartesian coordinate axis, n d Represents the dimension of the space, 2n d Indicates the total number of cube sampling points extracted;

[0151] Project the cube sampling point into the state space to obtain the cube sampling point X at time k-1 i,k-1∣k-1 :

[0152]

[0153] Among them, P k-1∣k-1 is the covariance matrix of the state estimate at time k-1, is the mean of the state estimate at time k-1;

[0154] Then X i,k-1∣k-1 Substituting into the nonlinear dynamic model, the cube after propagation adopts point That is, the predicted state at time k:

[0155]

[0156] Where f(·) is the discrete-time state prediction equation. The state error propagation equation is discretized by numerical integration to calculate the predicted state of the system from time k-1 to time k. The state error propagation equation of the system is:

[0157]

[0158] in, is the time rate of change of the pitch, roll, and heading attitude errors from time k-1 to time k, is the time rate of change of the velocity error from time k-1 to time k;

[0159] δΦ k-1 × indicates opposition to matrix form and has:

[0160] is the angular velocity of the navigation system at time k-1, is the angular velocity error in the navigation system at time k-1,

[0161] is the angular velocity of the Earth's rotation, is the error of the Earth's rotation angular velocity:

[0162]

[0163] Among them, Ω is the rotation rate of the earth, L is the latitude, and δL is the latitude error;

[0164] is the navigation error angular velocity at time k:

[0165]

[0166] Where R is the mean radius of curvature of the earth; h r is the height of the pipeline robot; are the eastward and northward velocities of the pipeline robot at time k-1 respectively; is the speed error of the pipeline robot in the east and north directions at time k-1; δh r is the height error of the pipeline robot;

[0167] N k-1 is the compensation term caused by the earth's rotation and navigation error at time k-1, M k-1 is the compensation term caused by the inertial system error at time k-1,

[0168] Calculate the predicted state mean:

[0169]

[0170] Compute the forecast error covariance:

[0171]

[0172] Process noise covariance matrix Q k-1 Adaptive adjustment based on IMU raw gyroscope data:

[0173] Improve the stability of the filter under different conditions, where Q0 is the initial value, and the empirical value 10 is taken in this invention. -3 , is the adaptive algorithm parameter, is the raw angular velocity measured by the gyroscope.

[0174] (24) Adaptive cubic Kalman filtering is used for measurement update, and the specific calculation steps are as follows:

[0175] (a) Calculate the pitch angle θ of the pipeline robot at time k by calculating the accelerometer and geomagnetic sensor data k , yaw angle ψ k , roll angle β k :

[0176]

[0177] Where μ is the relative magnetic permeability, They represent the acceleration components in the x, y, and z directions of the body coordinate system measured by the IMU accelerometer at time k, is the projection of the geomagnetic field at time k in the coordinate system of the pipeline robot, |B| is the total intensity of the geomagnetic field;

[0178] Calculate the transformation matrix from the body coordinate system to the navigation coordinate system at time k

[0179]

[0180] (b) Setting the threshold ξ g , determine whether the current IMU data is available:

[0181]

[0182] Among them, g z is the acceleration due to gravity, and the constant value of the present invention is 9.81 m / s 2 ,× represents the cross product,ξ g is a constant threshold, and the present invention takes an empirical value of 0.1.

[0183] (c) Obtain the observation vector based on the odometer, accelerometer, and magnetometer sensors: The speed measurement error is is the speed in the navigation coordinate system at time k, is the transformation matrix from the fuselage coordinate system to the navigation coordinate system at time k, is the speed in the fuselage coordinate system measured by the odometer at time k, is the projection of the geomagnetic field at time k in the coordinate system of the pipeline robot;

[0184] (d) Calculate the new cube sampling point X based on the predicted state covariance i,k∣k-1 :

[0185]

[0186] Among them, X i,k∣k-1 represents the i-th cubic point, represents the cubic point after the propagation of the predicted state covariance at time k-1, P k∣k-1 is the predicted state covariance matrix, which represents the uncertainty of the state variables, E{X k∣k-1} is the predicted state mean, which is used to adjust the center position of the cube point;

[0187] Calculate the cube sampling points P in the measurement space using the measurement model xz,k∣k-1 :

[0188] Z i,k∣k-1 =h(X i,k∣k-1) (13)

[0189] Where h(·) is the measurement model;

[0190] Calculates predicted measurement estimates

[0191]

[0192] Calculate the covariance matrix S for the uncertainty between the predicted and actual measurements k∣k-1 :

[0193]

[0194] Among them, R k is the measurement noise covariance matrix;

[0195] Calculate the state and measurement cross covariance matrix P xz,k∣k-1 :

[0196]

[0197] Calculate the Kalman gain K k :

[0198]

[0199] (25) Correct the pipeline robot status error:

[0200] (a) Update the state estimate and calculate the predicted state mean E{X k|k}:

[0201]

[0202] Among them, Z k is the actual measurement value of the sensor, is the predicted measurement value.

[0203] (b) Update the error covariance P k|k :

[0204]

[0205] Among them, P k∣k-1 is the predicted state covariance matrix, P zz,k∣k-1 Represents the forecast measurement error covariance matrix:

[0206]

[0207] (c) Establish the velocity error correction equation and attitude error correction equation, and update the system state parameters:

[0208] Establish the speed error correction equation:

[0209]

[0210] in, is the corrected velocity value at time k, is the speed value calculated at time k, [X k (5),X k (6),X k (7)] is the velocity error correction obtained by Kalman filtering at time k.

[0211] Establish the attitude error correction equation:

[0212]

[0213] in, is the modified k-time attitude transformation matrix, is the original k-1 moment posture transformation matrix, [X k (1),X k (2),X k (3)] ​​is the attitude error correction obtained by Kalman filtering. (26)

[0215] Update the position information of the pipeline robot:

[0216]

[0217] Among them, Pos k-1 represents the position information of the pipeline robot at time k-1, that is, the position at the previous moment, and Δt represents the time difference between time k and time k-1.

[0218] (3) Figure 3 As shown in the figure, an unsupervised pipeline anomaly detection method based on multi-scale feature fitting of a graphical model is used to analyze the image frame Image(k) collected by the camera to realize the detection of abnormal state inside the pipeline.

[0219] (35) Load the pre-trained ResNet-50 network model as a feature extractor;

[0220] (36) Construct the feature representation of the normal image set, that is, the image gallery set F G ;

[0221] (a) Load a set of M defect-free underground pipeline interior images as the normal image set D = {D1, D2, D3, ..., D M}, each image Among them, 1≤i≤M represents the total number of normal images, and w, h, c are the width, height, and length of the image respectively.

[0222] (b) ResNet-50 network is used to extract the feature map of the normal image set. Suppose the feature map extracted by the lth layer of the network is Among them, l∈{1,2,3,4}, w l ,h l Respectively represent the width and height of the feature map, c l Represents the number of channels (feature dimension) of the feature map of the lth layer. The feature vector of each pixel point p on the lth layer is expressed as 1≤p≤p l , where p l =w l ×h l Represents the total number of feature points in the feature map of layer l;

[0223] (c) For all normal images without defects D i ,1≤i≤M and all layers l∈{1,2,3,4}, summarize the feature vectors of all feature points to form a local gallery set

[0224]

[0225] (d) Perform global average pooling on each image to obtain the global gallery set

[0226]

[0227] (e) Combine the local gallery set with the global gallery set to construct the final complete gallery set

[0228] (33) Use the ResNet-50 network to collect images in real time from the pipeline robot camera Perform feature extraction, w, h, c are the width, height and length of the image respectively, and construct the query set F Q ;

[0229] (a) Using ResNet-50 network to extract real-time acquired images The feature map of the network is as follows: Among them, l∈{1,2,3,4}, w l ,h l Respectively represent the width and height of the feature map, c l Represents the number of channels (feature dimension) of the feature map of the lth layer. The feature vector of each pixel point p on the lth layer is expressed as 1≤p≤p l , where p l =wl ×h l Represents the total number of feature points in the feature map of layer l.

[0230] (b) For Image(k) and all layers l∈{1,2,3,4}, the feature vectors of all feature points are aggregated to form a local query set

[0231]

[0232] (c) Perform global average pooling on each image to obtain the global query set

[0233]

[0234] (d) Combine the local query set with the global query set to construct the final complete query set

[0235] (34) For each query point y i ∈F Q , in Gallery Set F G Find k s Nearest neighbors {f i}, and calculate the cosine similarity a i,j :

[0236]

[0237] Select the first K (in this invention, the value is K = 5) f with the maximum similarity j As y i The nearest neighbor, that is, N i = {v g,1 ,v g,2 ,…,v g,j ,…,v g,K}, the vertices corresponding to these nearest neighbor features are represented as v g,j ∈V G , V G is a set of image gallery vertices (representing local feature points extracted from normal training images), and constructs an adjacency matrix A, with a i,j ∈A:

[0238]

[0239] (35) For query point y i , using K nearest neighbors for weighted feature fitting:

[0240]

[0241] Calculate the anomaly score s for each vertex based on the cosine distancei :

[0242]

[0243] Get the anomaly score matrix And use bilinear interpolation method to Zoom in to the same size as the original test image:

[0244]

[0245] According to the weights of different layers, the weighted multi-scale anomaly score matrix Ω(Image(k)) is calculated:

[0246]

[0247] in, Indicates the weighting coefficient of the feature map of the first layer. In the present invention, the empirical coefficients are: w1=0.1, w2=0.2, w3=0.3, w4=0.4.

[0248] (38) Set a threshold value τ1 (in the present invention, τ1 takes the empirical value of 0.85) to determine whether the pixel is abnormal. If Ω(Image(k),p)>τ1, the pixel point p is marked as an abnormal pixel. Set a threshold value τ2 (in the present invention, τ2 takes the empirical value of 0.95). If the maximum abnormal value of the entire image max(Ω(Image(k)))>τ2, Image(k) is judged to be an abnormal image, and a binary abnormal segmentation map Defect(Image(k)) is generated:

[0249]

[0250] (4) When an abnormal situation is detected in the pipeline through step (3), the current robot positioning coordinate Pos is fed back k , to achieve accurate feedback of abnormal position {Image(k),Defect(Image(k)),Pos k}.

Claims

1. A method for locating abnormalities inside underground pipelines by combining inertial navigation and vision, characterized in that: The method comprises the following steps: (1) An active crawling pipeline robot is equipped with a camera, a fill light, an inertial measurement unit (IMU) with a built-in three-axis accelerometer and three-axis gyroscope, a wheel odometer, and a geomagnetic sensor to collect and transmit data inside the pipeline; (2) The positioning method based on adaptive cubic Kalman filter is used to perform real-time positioning calculation of the pipeline robot and provide the coordinate data Pos of the pipeline robot at each moment. k ; (3) An unsupervised pipeline anomaly detection method based on multi-scale feature fitting of a graphical model is used to analyze the image frame Image(k) collected by the camera to detect abnormal conditions inside the pipeline. (4) When an abnormal situation is detected in the pipeline through step (3), the current robot positioning coordinate Pos is fed back k , achieving accurate feedback of abnormal positions.

2. The method for locating anomalies inside underground pipelines by combined inertial navigation and vision collaboration according to claim 1 is characterized in that: The step (1) specifically comprises the following steps: (11) Before the pipeline robot is put into operation, the sensor needs to be calibrated. First, the accelerometer and gyroscope of the inertial measurement unit (IMU) are bias-corrected, and the initial orientation of the geomagnetic sensor is determined using static measurement. (12) After the pipeline robot enters the pipeline, it turns on the camera and fill light, and starts to collect the robot's own positioning data and image data. Assume that the data collected at time k are: The acceleration measurement value of the inertial measurement unit IMU in the body coordinate system is They represent the acceleration measurement values ​​of the pipeline robot in the x-axis, y-axis, and z-axis directions in the body coordinate system at time k, respectively. The superscript T represents the transpose of the vector and the angular velocity measurement value of the inertial measurement unit IMU in the body coordinate system. They represent the angular velocity measurement values ​​of the pipeline robot in the x-axis, y-axis, and z-axis directions in the body coordinate system at time k, and the odometer measurement speed in the body coordinate system δs k Indicates the speed measured by the odometer in the forward or backward direction of the pipeline robot and the geomagnetic measurement value in the body coordinate system They represent the geomagnetic measurement values ​​of the pipeline robot at the east, north, and up directions at time k, and the data obtained by the gravity accelerometer g = [g x ,g y ,g z ] T , the image frame Image(k) acquired by the camera; where the subscripts x, y, and z respectively represent the "x, y, and z" axes in the body coordinate system, pointing to the east, north, and up directions respectively; the subscript b represents the body coordinate system of the pipeline robot; it is calculated by measuring the rotation speed of the wheels; "k" represents the kth moment, and when k = 0, it represents the initial state value; (13) The collected data will be transmitted to the ground workstation for processing via the optical fiber cable at the rear end of the active crawling pipeline robot.

3. The method for locating anomalies inside underground pipelines by combined inertial navigation and vision collaboration according to claim 2 is characterized in that: The step (2) specifically comprises the following steps: (21) Define the state vector of the pipeline robot at time k as X k , Definition Φ k is the pipeline robot pose, are the pitch, roll and yaw angles of the pipeline robot at time k; δΦ k is the navigation angle error of the robot at time k, in, They represent the pitch, roll and yaw angle errors of the pipeline robot at time k respectively; V k represents the speed of the pipeline robot at time k in the navigation coordinate system, They represent the speed of the pipeline robot in the east, north and upward directions at time k; δV k is the velocity error of the pipeline robot in the navigation coordinate system at time k, They represent the velocity errors of the pipeline robot in the east, north, and up directions at time k, respectively; is the bias error of the angular velocity measurement value of the inertial measurement unit IMU of the pipeline robot in the body coordinate system at time k, They respectively represent the bias errors of the angular velocity measurement values ​​of the inertial measurement unit IMU in the x-axis, y-axis, and z-axis directions in the body coordinate system of the pipeline robot at time k; represents the IMU accelerometer bias error of the pipeline robot in the body coordinate system at time k, They represent the IMU accelerometer bias errors in the x-axis, y-axis, and z-axis directions of the pipeline robot in the body coordinate system at time k, respectively; Define the state error covariance matrix of the pipeline robot at time k as P k , They respectively represent the variance of the pipeline robot's posture at time k, the velocity variance in the navigation coordinate system, the bias error variance of the angular velocity measurement value of the inertial measurement unit IMU of the pipeline robot in the body coordinate system at time k, and the bias error variance of the IMU accelerometer in the body coordinate system; Define the process noise covariance matrix of the pipeline robot at time k as Q k , They respectively represent the variance of the angular velocity measurement value of the inertial measurement unit IMU in the body coordinate system and the variance of the acceleration measurement value of the inertial measurement unit IMU in the body coordinate system; Define the measurement noise covariance matrix of the pipeline robot at time k as R k , R k = They respectively represent the variance of the speed measured by the odometer in the body coordinate system, the variance of the geomagnetic measurement value in the body coordinate system, and the variance of the acceleration measurement value of the inertial measurement unit IMU in the body coordinate system; Define the coordinates of the pipeline robot at time k as Pos k ,Pos k =[x k ,y k ,z k ] m , x k ,y k 、z k They represent the coordinates of the x, y, and z axes in the body coordinate system, respectively, and the superscript T represents the transpose of the vector; definition is the transformation matrix from the body coordinate system to the navigation coordinate system at time k; (22) Pipeline robot sensor data initialization: Initialize state vector: Initialize the state error covariance matrix: Initialize the process noise covariance matrix: Initialize the measurement noise covariance matrix: Initialize the pipeline robot positioning coordinates: Pos0 = [x0, y0, z0] T , the transformation matrix of the initial posture from the body coordinate system to the navigation coordinate system is defined as: (23) Adaptive cubic Kalman filter is used for state prediction, and the specific calculation steps are as follows: From d Draw cube sampling points at the intersection of the dimensional unit sphere and the Cartesian coordinate axes and perform a scale transformation: Among them, e i is a unit vector, the i-1th element is 1, the rest are 0, ξ i Indicates that from n d The result of scaling after the cube sampling points are extracted from the intersection of the unit sphere and the Cartesian coordinate axis, n d Represents the dimension of the space, 2n d Indicates the total number of cube sampling points extracted; Project the cube sampling point into the state space to obtain the cube sampling point X at time k-1 i,k-1∣k-1 : Among them, P k-1∣k-1 is the covariance matrix of the state estimate at time k-1, is the mean of the state estimate at time k-1; Then X i,k-1∣k-1 Substituting into the nonlinear dynamic model, the cube after propagation adopts point That is, the predicted state at time k: Where f(·) is the discrete-time state prediction equation. The state error propagation equation is discretized by numerical integration to calculate the predicted state of the system from time k-1 to time k. The state error propagation equation of the system is: in, is the time rate of change of the pitch, roll, and heading attitude errors from time k-1 to time k, is the time rate of change of the velocity error from time k-1 to time k; δΦ k-1 × indicates opposition to matrix form and has: is the angular velocity of the navigation system at time k-1, is the angular velocity error in the navigation system at time k-1, is the angular velocity of the Earth's rotation, is the error of the Earth's rotation angular velocity: Among them, Ω is the rotation rate of the earth, L is the latitude, and δL is the latitude error; is the navigation error angular velocity at time k: Where R is the mean radius of curvature of the earth; h r is the height of the pipeline robot; are the eastward and northward velocities of the pipeline robot at time k-1 respectively; and is the speed error of the pipeline robot in the east and north directions at time k-1; δh r is the height error of the pipeline robot; N k-1 is the compensation term caused by the earth's rotation and navigation error at time k-1, M k-1 is the compensation term caused by the inertial system error at time k-1, Calculate the predicted state mean: Compute the forecast error covariance: Process noise covariance matrix Q k-1 Adaptive adjustment based on IMU raw gyroscope data: Improve the stability of the filter in different states, where Q0 is the initial value, is the adaptive algorithm parameter, is the original angular velocity measured by the gyroscope; (24) Adaptive cubic Kalman filtering is used for measurement update, and the specific calculation steps are as follows: (a) Calculate the pitch angle θ of the pipeline robot at time k by calculating the accelerometer and geomagnetic sensor data k , yaw angle ψ k , roll angle β k : Where μ is the relative magnetic permeability, They represent the acceleration components in the x, y, and z directions in the body coordinate system measured by the IMU accelerometer at time k, is the projection of the geomagnetic field at time k in the coordinate system of the pipeline robot, |B| is the total intensity of the geomagnetic field; Calculate the transformation matrix from the body coordinate system to the navigation coordinate system at time k (b) Setting the threshold ξ g , determine whether the current IMU data is available: Among them, g z is the acceleration due to gravity, × represents the cross product, ξ g is a constant threshold; (c) Obtain the observation vector based on the odometer, accelerometer, and magnetometer sensors: The speed measurement error is is the speed in the navigation coordinate system at time k, is the transformation matrix from the fuselage coordinate system to the navigation coordinate system at time k, is the speed in the fuselage coordinate system measured by the odometer at time k, is the projection of the geomagnetic field at time k in the coordinate system of the pipeline robot; (d) Calculate the new cube sampling point X based on the predicted state covariance i,k|k-1 : Among them, X i，k|k-1 represents the i-th cubic point, represents the cubic point after the propagation of the predicted state covariance at time k-1, P k|k-1 is the predicted state covariance matrix, which represents the uncertainty of the state variables, E{X k|k-1 } is the predicted state mean, which is used to adjust the center position of the cube point; Calculate the cube sampling points P in the measurement space using the measurement model xz，k|k-1 : Z i，k|k-1 =h(X i,k|k-1 ) (13) Where h(·) is the measurement model; Calculates predicted measurement estimates Calculate the covariance matrix S for the uncertainty between the predicted and actual measurements k|k-1 : Among them, R k is the measurement noise covariance matrix; Calculate the state and measurement cross covariance matrix P xz，k|k-1 : Calculate the Kalman gain K k : (25) Correct the pipeline robot status error: (a) Update the state estimate and calculate the predicted state mean E{X k|k }: Among them, Z k is the actual measurement value of the sensor, is the predicted measurement; (b) Update the error covariance P k|k : P k|k =P k|k-1 -K k P zz，k|k-1 (K k ) T (19) Among them, P k|k-1 is the predicted state covariance matrix, P zz，k|k-1 Represents the forecast measurement error covariance matrix: (c) Establish the velocity error correction equation and attitude error correction equation, and update the system state parameters: Establish the speed error correction equation: in, is the corrected velocity value at time k, is the speed value calculated at time k, [X k (5), X k (6), X k (7)] is the velocity error correction obtained by Kalman filtering at time k; Establish the attitude error correction equation: in, is the modified k-time attitude transformation matrix, is the original k-1 moment posture transformation matrix, [X k (1), X k (2), X k (3)] ​​is the attitude error correction obtained by Kalman filtering; (26) Update the location information of the pipeline robot: Among them, Pos k-1 represents the position information of the pipeline robot at time k-1, that is, the position at the previous moment, and Δt represents the time difference between time k and time k-1.

4. The method for locating anomalies inside underground pipelines by combined inertial navigation and vision collaboration according to claim 3 is characterized in that: The step (3) specifically comprises the following steps: (31) Load the pre-trained ResNet-50 network model as a feature extractor; (32) Construct the feature representation of the normal image set, that is, the image gallery set F G ; (a) Load a set of M defect-free underground pipeline interior images as the normal image set D = {D1, D2, D3, ..., D M }, each image Among them, 1≤i≤M represents the total number of normal images, w, h, c are the width, height and length of the image respectively; (b) ResNet-50 network is used to extract the feature map of the normal image set. Suppose the feature map extracted by the lth layer of the network is Among them, l∈{1, 2, 3, 4}, w l ,h l Respectively represent the width and height of the feature map, c l Represents the number of channels of the feature map of the lth layer. The feature vector of each pixel point p on the lth layer is expressed as where p l =w l ×h l Represents the total number of feature points in the feature map of layer l; (c) For all normal images without defects D i , 1≤i≤M and all layers l∈{1, 2, 3, 4}, summarize the feature vectors of all feature points to form a local gallery set (d) Perform global average pooling on each image to obtain the global gallery set (e) Combine the local gallery set with the global gallery set to construct the final complete gallery set (33) Use the ResNet-50 network to collect images in real time from the pipeline robot camera Perform feature extraction, w, h, c are the width, height and length of the image respectively, and construct the query set F Q ; (a) Using ResNet-50 network to extract real-time acquired images The feature map of the network is as follows: Among them, l∈{1, 2, 3, 4}, w i ,h l Respectively represent the width and height of the feature map, c l Represents the number of channels of the feature map of the lth layer; the feature vector of each pixel point p on the lth layer is expressed as where p l =w l ×h l Represents the total number of feature points in the feature map of layer l; (b) For Image(k) and all layers l∈{1, 2, 3, 4}, summarize the feature vectors of all feature points to form a local query set (c) Perform global average pooling on each image to obtain the global query set (d) Combine the local query set with the global query set to construct the final complete query set (34) For each query point y i ∈F Q , in Gallery Set F G Find k s Nearest neighbors {f i }, and calculate the cosine similarity a i，j : Select the first K largest similarity f j As y i The nearest neighbor, that is, N i = {v g,1 , v g,2 , …, v g，j , …, v g,K }, the vertices corresponding to these nearest neighbor features are represented as v g，j ∈V G , V G is a set of graph vertices, and constructs an adjacency matrix A, with a i，j ∈A: (35) For query point y i , using K nearest neighbors for weighted feature fitting: Calculate the anomaly score s for each vertex based on the cosine distance i : Get the anomaly score matrix And use bilinear interpolation method to Zoom in to the same size as the original test image: According to the weights of different layers, the weighted multi-scale anomaly score matrix Ω(Image(k)) is calculated: in, Represents the weight coefficient of the feature map of the lth layer; (36) Set the threshold τ1 to determine whether the pixel is abnormal. If Ω(Image(k), p)>τ1, the pixel p is marked as an abnormal pixel. Set the threshold τ2. If the maximum abnormal value of the entire image max(Ω(Image(k)))>τ2, Image(k) is judged as an abnormal image, and a binary abnormal segmentation map Defect(Image(k)) is generated:

Images