IT评测·应用市场-qidao123.com技术社区

标题: LIO-SAM代码解析:imuPreintegration.cpp(一) [打印本页]

作者: 水军大提督    时间: 2025-1-12 10:09
标题: LIO-SAM代码解析:imuPreintegration.cpp(一)

1.imuHandler

imuHandler 函数用于处理IMU(惯性丈量单元)数据,通过积分IMU丈量值预测当前的位姿与速度,并将IMU位姿转换到激光雷达坐标系,最终发布惯导里程计数据。
  1. // IMU 数据处理回调函数
  2. void imuHandler(const sensor_msgs::Imu::ConstPtr& imu_raw)
  3. {
  4.     // 加锁,确保线程安全
  5.     std::lock_guard<std::mutex> lock(mtx);
  6.     // 转换原始 IMU 数据
  7.     sensor_msgs::Imu thisImu = imuConverter(*imu_raw);
  8.     // 将转换后的 IMU 数据加入优化队列和 IMU 队列
  9.     imuQueOpt.push_back(thisImu);
  10.     imuQueImu.push_back(thisImu);
  11.     // 如果尚未完成第一次优化,则退出函数
  12.     if (doneFirstOpt == false)
  13.         return;
  14.     // 获取当前 IMU 数据的时间戳
  15.     double imuTime = ROS_TIME(&thisImu);
  16.    
  17.     // 计算时间间隔 dt
  18.     double dt = (lastImuT_imu < 0) ? (1.0 / 500.0) : (imuTime - lastImuT_imu);
  19.     lastImuT_imu = imuTime;
  20.     // 积分当前 IMU 测量数据
  21.     imuIntegratorImu_->integrateMeasurement(
  22.         gtsam::Vector3(thisImu.linear_acceleration.x, thisImu.linear_acceleration.y, thisImu.linear_acceleration.z),
  23.         gtsam::Vector3(thisImu.angular_velocity.x,    thisImu.angular_velocity.y,    thisImu.angular_velocity.z),
  24.         dt
  25.     );
  26.     // 使用积分后的状态预测当前的导航状态
  27.     gtsam::NavState currentState = imuIntegratorImu_->predict(prevStateOdom, prevBiasOdom);
  28.     // 创建并初始化里程计消息
  29.     nav_msgs::Odometry odometry;
  30.     odometry.header.stamp = thisImu.header.stamp;
  31.     odometry.header.frame_id = odometryFrame;
  32.     odometry.child_frame_id = "odom_imu";
  33.     // 将 IMU 位姿转换到激光雷达坐标系
  34.     gtsam::Pose3 imuPose = gtsam::Pose3(currentState.quaternion(), currentState.position());
  35.     gtsam::Pose3 lidarPose = imuPose.compose(imu2Lidar);
  36.     // 填充里程计位置信息
  37.     odometry.pose.pose.position.x = lidarPose.translation().x();
  38.     odometry.pose.pose.position.y = lidarPose.translation().y();
  39.     odometry.pose.pose.position.z = lidarPose.translation().z();
  40.     odometry.pose.pose.orientation.x = lidarPose.rotation().toQuaternion().x();
  41.     odometry.pose.pose.orientation.y = lidarPose.rotation().toQuaternion().y();
  42.     odometry.pose.pose.orientation.z = lidarPose.rotation().toQuaternion().z();
  43.     odometry.pose.pose.orientation.w = lidarPose.rotation().toQuaternion().w();
  44.    
  45.     // 填充里程计速度信息
  46.     odometry.twist.twist.linear.x = currentState.velocity().x();
  47.     odometry.twist.twist.linear.y = currentState.velocity().y();
  48.     odometry.twist.twist.linear.z = currentState.velocity().z();
  49.     odometry.twist.twist.angular.x = thisImu.angular_velocity.x + prevBiasOdom.gyroscope().x();
  50.     odometry.twist.twist.angular.y = thisImu.angular_velocity.y + prevBiasOdom.gyroscope().y();
  51.     odometry.twist.twist.angular.z = thisImu.angular_velocity.z + prevBiasOdom.gyroscope().z();
  52.    
  53.     // 发布里程计消息
  54.     pubImuOdometry.publish(odometry);
  55. }
复制代码
1. IMU丈量值积分
IMU数据包含线加速度                                              a                            imu                                       \mathbf{a}_{\text{imu}}                  aimu​ 和角速度                                              ω                            imu                                       \boldsymbol{\omega}_{\text{imu}}                  ωimu​,函数通过GTSAM库的 integrateMeasurement 方法对IMU丈量值进行积分,公式为:
                                                                                           v                                                       k                                           +                                           1                                                                                                                                     =                                                       v                                           k                                                      +                                        (                                                       R                                           k                                                      ⋅                                        (                                                       a                                           imu                                                      −                                                       b                                           a                                                      )                                        +                                        g                                        )                                        ⋅                                        Δ                                        t                                        ,                                                                                                                            p                                                       k                                           +                                           1                                                                                                                                     =                                                       p                                           k                                                      +                                                       v                                           k                                                      ⋅                                        Δ                                        t                                        +                                                       1                                           2                                                      (                                                       R                                           k                                                      ⋅                                        (                                                       a                                           imu                                                      −                                                       b                                           a                                                      )                                        +                                        g                                        )                                        ⋅                                        Δ                                                       t                                           2                                                      ,                                                                                                                            R                                                       k                                           +                                           1                                                                                                                                     =                                                       R                                           k                                                      ⋅                                        exp                                        ⁡                                        (                                        (                                                       ω                                           imu                                                      −                                                       b                                           ω                                                      )                                        ⋅                                        Δ                                        t                                        )                                        ,                                                                                \begin{aligned} \mathbf{v}_{k+1} &= \mathbf{v}_k + (\mathbf{R}_k \cdot (\mathbf{a}_{\text{imu}} - \mathbf{b}_a) + \mathbf{g}) \cdot \Delta t, \\ \mathbf{p}_{k+1} &= \mathbf{p}_k + \mathbf{v}_k \cdot \Delta t + \frac{1}{2} (\mathbf{R}_k \cdot (\mathbf{a}_{\text{imu}} - \mathbf{b}_a) + \mathbf{g}) \cdot \Delta t^2, \\ \mathbf{R}_{k+1} &= \mathbf{R}_k \cdot \exp((\boldsymbol{\omega}_{\text{imu}} - \mathbf{b}_\omega) \cdot \Delta t), \end{aligned}                     vk+1​pk+1​Rk+1​​=vk​+(Rk​⋅(aimu​−ba​)+g)⋅Δt,=pk​+vk​⋅Δt+21​(Rk​⋅(aimu​−ba​)+g)⋅Δt2,=Rk​⋅exp((ωimu​−bω​)⋅Δt),​
此中:


2. 位姿预测
通过上述积分盘算,预测得到当前状态:
                                                    State                                           k                                  +                                  1                                                 =                                       [                                                                                                     p                                                               k                                                 +                                                 1                                                                                                                                                                  v                                                               k                                                 +                                                 1                                                                                                                                                                  R                                                               k                                                 +                                                 1                                                                                                                ]                                      .                                  \text{State}_{k+1} = \begin{bmatrix} \mathbf{p}_{k+1} \\ \mathbf{v}_{k+1} \\ \mathbf{R}_{k+1} \end{bmatrix}.                     Statek+1​=               ​pk+1​vk+1​Rk+1​​               ​.
函数中调用 imuIntegratorImu_->predict 实现此预测。

3. 位姿转换(IMU到激光雷达坐标系)
通过以下公式将IMU的位姿                                              T                            imu                                       \mathbf{T}_{\text{imu}}                  Timu​ 转换为激光雷达位姿                                              T                            lidar                                       \mathbf{T}_{\text{lidar}}                  Tlidar​:
                                                    T                               lidar                                      =                                       T                               imu                                      ⋅                                       T                               imu2lidar                                      ,                                  \mathbf{T}_{\text{lidar}} = \mathbf{T}_{\text{imu}} \cdot \mathbf{T}_{\text{imu2lidar}},                     Tlidar​=Timu​⋅Timu2lidar​,
此中                                              T                            imu2lidar                                       \mathbf{T}_{\text{imu2lidar}}                  Timu2lidar​ 是IMU到激光雷达的静态外参。
对于位姿                                    T                              \mathbf{T}                  T,其包含旋转                                    R                              \mathbf{R}                  R 宁静移                                    t                              \mathbf{t}                  t,即:
                                         T                            =                                       [                                                                                     R                                                                                             t                                                                                                                   0                                                                                             1                                                                                 ]                                      .                                  \mathbf{T} = \begin{bmatrix} \mathbf{R} & \mathbf{t} \\ \mathbf{0} & 1 \end{bmatrix}.                     T=[R0​t1​].

4. 惯导里程计发布
最终通过以下公式将预测得到的激光雷达位姿和速度发布为惯导里程计数据:


2. odometryHandler

odometryHandler 函数的重要功能是接收激光里程计数据,结合IMU预积分对系统状态进行优化,并更新系统的位姿、速度和偏置。以下为具体的数学过程。
  1. void odometryHandler(const nav_msgs::Odometry::ConstPtr& odomMsg)
  2. {
  3.     std::lock_guard<std::mutex> lock(mtx);  // 保护共享资源,防止多线程访问时出现竞争条件
  4.     double currentCorrectionTime = ROS_TIME(odomMsg);  // 获取当前里程计消息的时间戳
  5.     // 确保有IMU数据可供积分
  6.     if (imuQueOpt.empty())
  7.         return;
  8.     // 从里程计消息中提取位置和姿态信息
  9.     float p_x = odomMsg->pose.pose.position.x;
  10.     float p_y = odomMsg->pose.pose.position.y;
  11.     float p_z = odomMsg->pose.pose.position.z;
  12.     float r_x = odomMsg->pose.pose.orientation.x;
  13.     float r_y = odomMsg->pose.pose.orientation.y;
  14.     float r_z = odomMsg->pose.pose.orientation.z;
  15.     float r_w = odomMsg->pose.pose.orientation.w;
  16.     bool degenerate = (int)odomMsg->pose.covariance[0] == 1 ? true : false;  // 检查是否为退化情况
  17.     gtsam::Pose3 lidarPose = gtsam::Pose3(gtsam::Rot3::Quaternion(r_w, r_x, r_y, r_z), gtsam::Point3(p_x, p_y, p_z));  // 构造激光雷达的位姿
  18.     // 0. 初始化系统
  19.     if (systemInitialized == false)
  20.     {
  21.         resetOptimization();  // 重置优化器
  22.         // 弹出旧的IMU消息
  23.         while (!imuQueOpt.empty())
  24.         {
  25.             if (ROS_TIME(&imuQueOpt.front()) < currentCorrectionTime - delta_t)
  26.             {
  27.                 lastImuT_opt = ROS_TIME(&imuQueOpt.front());
  28.                 imuQueOpt.pop_front();
  29.             }
  30.             else
  31.                 break;
  32.         }
  33.         // 初始位姿
  34.         prevPose_ = lidarPose.compose(lidar2Imu);  // 将激光雷达位姿转换为IMU位姿
  35.         gtsam::PriorFactor<gtsam::Pose3> priorPose(X(0), prevPose_, priorPoseNoise);  // 添加初始位姿先验因子
  36.         graphFactors.add(priorPose);
  37.         // 初始速度
  38.         prevVel_ = gtsam::Vector3(0, 0, 0);  // 初始速度为零
  39.         gtsam::PriorFactor<gtsam::Vector3> priorVel(V(0), prevVel_, priorVelNoise);  // 添加初始速度先验因子
  40.         graphFactors.add(priorVel);
  41.         // 初始偏差
  42.         prevBias_ = gtsam::imuBias::ConstantBias();  // 初始偏差为零
  43.         gtsam::PriorFactor<gtsam::imuBias::ConstantBias> priorBias(B(0), prevBias_, priorBiasNoise);  // 添加初始偏差先验因子
  44.         graphFactors.add(priorBias);
  45.         // 添加初始值
  46.         graphValues.insert(X(0), prevPose_);
  47.         graphValues.insert(V(0), prevVel_);
  48.         graphValues.insert(B(0), prevBias_);
  49.         // 优化一次
  50.         optimizer.update(graphFactors, graphValues);
  51.         graphFactors.resize(0);
  52.         graphValues.clear();
  53.         imuIntegratorImu_->resetIntegrationAndSetBias(prevBias_);  // 重置IMU积分器并设置偏差
  54.         imuIntegratorOpt_->resetIntegrationAndSetBias(prevBias_);
  55.         key = 1;
  56.         systemInitialized = true;  // 标记系统已初始化
  57.         return;
  58.     }
  59.     // 重置图以提高速度
  60.     if (key == 100)
  61.     {
  62.         // 获取重置前的噪声模型
  63.         gtsam::noiseModel::Gaussian::shared_ptr updatedPoseNoise = gtsam::noiseModel::Gaussian::Covariance(optimizer.marginalCovariance(X(key-1)));
  64.         gtsam::noiseModel::Gaussian::shared_ptr updatedVelNoise  = gtsam::noiseModel::Gaussian::Covariance(optimizer.marginalCovariance(V(key-1)));
  65.         gtsam::noiseModel::Gaussian::shared_ptr updatedBiasNoise = gtsam::noiseModel::Gaussian::Covariance(optimizer.marginalCovariance(B(key-1)));
  66.         // 重置图
  67.         resetOptimization();
  68.         // 添加位姿
  69.         gtsam::PriorFactor<gtsam::Pose3> priorPose(X(0), prevPose_, updatedPoseNoise);
  70.         graphFactors.add(priorPose);
  71.         // 添加速度
  72.         gtsam::PriorFactor<gtsam::Vector3> priorVel(V(0), prevVel_, updatedVelNoise);
  73.         graphFactors.add(priorVel);
  74.         // 添加偏差
  75.         gtsam::PriorFactor<gtsam::imuBias::ConstantBias> priorBias(B(0), prevBias_, updatedBiasNoise);
  76.         graphFactors.add(priorBias);
  77.         // 添加值
  78.         graphValues.insert(X(0), prevPose_);
  79.         graphValues.insert(V(0), prevVel_);
  80.         graphValues.insert(B(0), prevBias_);
  81.         // 优化一次
  82.         optimizer.update(graphFactors, graphValues);
  83.         graphFactors.resize(0);
  84.         graphValues.clear();
  85.         key = 1;
  86.     }
  87.     // 1. 积分IMU数据并优化
  88.     while (!imuQueOpt.empty())
  89.     {
  90.         // 弹出并积分两个优化之间的IMU数据
  91.         sensor_msgs::Imu *thisImu = &imuQueOpt.front();
  92.         double imuTime = ROS_TIME(thisImu);
  93.         if (imuTime < currentCorrectionTime - delta_t)
  94.         {
  95.             double dt = (lastImuT_opt < 0) ? (1.0 / 500.0) : (imuTime - lastImuT_opt);
  96.             imuIntegratorOpt_->integrateMeasurement(
  97.                     gtsam::Vector3(thisImu->linear_acceleration.x, thisImu->linear_acceleration.y, thisImu->linear_acceleration.z),
  98.                     gtsam::Vector3(thisImu->angular_velocity.x,    thisImu->angular_velocity.y,    thisImu->angular_velocity.z), dt);
  99.             
  100.             lastImuT_opt = imuTime;
  101.             imuQueOpt.pop_front();
  102.         }
  103.         else
  104.             break;
  105.     }
  106.     // 向图中添加IMU因子
  107.     const gtsam::PreintegratedImuMeasurements& preint_imu = dynamic_cast<const gtsam::PreintegratedImuMeasurements&>(*imuIntegratorOpt_);
  108.     gtsam::ImuFactor imu_factor(X(key - 1), V(key - 1), X(key), V(key), B(key - 1), preint_imu);
  109.     graphFactors.add(imu_factor);
  110.     // 添加IMU偏差之间的因子
  111.     graphFactors.add(gtsam::BetweenFactor<gtsam::imuBias::ConstantBias>(B(key - 1), B(key), gtsam::imuBias::ConstantBias(),
  112.                          gtsam::noiseModel::Diagonal::Sigmas(sqrt(imuIntegratorOpt_->deltaTij()) * noiseModelBetweenBias)));
  113.     // 添加位姿因子
  114.     gtsam::Pose3 curPose = lidarPose.compose(lidar2Imu);
  115.     gtsam::PriorFactor<gtsam::Pose3> pose_factor(X(key), curPose, degenerate ? correctionNoise2 : correctionNoise);
  116.     graphFactors.add(pose_factor);
  117.     // 插入预测值
  118.     gtsam::NavState propState_ = imuIntegratorOpt_->predict(prevState_, prevBias_);
  119.     graphValues.insert(X(key), propState_.pose());
  120.     graphValues.insert(V(key), propState_.v());
  121.     graphValues.insert(B(key), prevBias_);
  122.     // 优化
  123.     optimizer.update(graphFactors, graphValues);
  124.     optimizer.update();
  125.     graphFactors.resize(0);
  126.     graphValues.clear();
  127.     // 为下一步重置预积分
  128.     gtsam::Values result = optimizer.calculateEstimate();
  129.     prevPose_  = result.at<gtsam::Pose3>(X(key));
  130.     prevVel_   = result.at<gtsam::Vector3>(V(key));
  131.     prevState_ = gtsam::NavState(prevPose_, prevVel_);
  132.     prevBias_  = result.at<gtsam::imuBias::ConstantBias>(B(key));
  133.     // 重置优化预积分对象
  134.     imuIntegratorOpt_->resetIntegrationAndSetBias(prevBias_);
  135.     // 检查优化是否失败
  136.     if (failureDetection(prevVel_, prevBias_))
  137.     {
  138.         resetParams();
  139.         return;
  140.     }
  141.     // 2. 优化后,重新传播IMU里程计预积分
  142.     prevStateOdom = prevState_;
  143.     prevBiasOdom  = prevBias_;
  144.         // 首先弹出比当前校正数据旧的IMU消息
  145.         double lastImuQT = -1;
  146.         while (!imuQueImu.empty() && ROS_TIME(&imuQueImu.front()) < currentCorrectionTime - delta_t)
  147.         {
  148.             lastImuQT = ROS_TIME(&imuQueImu.front());
  149.             imuQueImu.pop_front();
  150.         }
  151.         // 重新传播
  152.         if (!imuQueImu.empty())
  153.         {
  154.             // 使用新优化的偏差重置偏差
  155.             imuIntegratorImu_->resetIntegrationAndSetBias(prevBiasOdom);
  156.             // 从本次优化的开始处积分IMU消息
  157.             for (int i = 0; i < (int)imuQueImu.size(); ++i)
  158.             {
  159.                 sensor_msgs::Imu *thisImu = &imuQueImu[i];
  160.                 double imuTime = ROS_TIME(thisImu);
  161.                 double dt = (lastImuQT < 0) ? (1.0 / 500.0) : (imuTime - lastImuQT);
  162.                 imuIntegratorImu_->integrateMeasurement(gtsam::Vector3(thisImu->linear_acceleration.x, thisImu->linear_acceleration.y, thisImu->linear_acceleration.z),
  163.                                                         gtsam::Vector3(thisImu->angular_velocity.x,    thisImu->angular_velocity.y,    thisImu->angular_velocity.z), dt);
  164.                 lastImuQT = imuTime;
  165.             }
  166.         }
  167.         ++key;
  168.         doneFirstOpt = true;
  169.     }
复制代码
1. 里程计与IMU的初始对齐
接收到的里程计位姿由以下参数表现:

激光雷达的位姿用 GTSAM 表现为:
                                                    T                               lidar                                      =                                       [                                                                                                     R                                              lidar                                                                                                                            p                                              lidar                                                                                                                                  0                                                                                             1                                                                                 ]                                      ,                                  \mathbf{T}_{\text{lidar}} = \begin{bmatrix} \mathbf{R}_{\text{lidar}} & \mathbf{p}_{\text{lidar}} \\ \mathbf{0} & 1 \end{bmatrix},                     Tlidar​=[Rlidar​0​plidar​1​],
此中                                              R                            lidar                                       \mathbf{R}_{\text{lidar}}                  Rlidar​ 由四元数                                              q                            lidar                                       \mathbf{q}_{\text{lidar}}                  qlidar​ 转换得到。
激光雷达位姿转换为IMU位姿:
                                                    T                               imu                                      =                                       T                               lidar                                      ⋅                                       T                               lidar2imu                                      ,                                  \mathbf{T}_{\text{imu}} = \mathbf{T}_{\text{lidar}} \cdot \mathbf{T}_{\text{lidar2imu}},                     Timu​=Tlidar​⋅Tlidar2imu​,
此中                                              T                            lidar2imu                                       \mathbf{T}_{\text{lidar2imu}}                  Tlidar2imu​ 为外参。
2. 系统初始化
初始位姿、速度和偏置,初始位姿                                              T                            0                                       \mathbf{T}_0                  T0​:
                                                    T                               0                                      =                                       T                               imu                                      .                                  \mathbf{T}_0 = \mathbf{T}_{\text{imu}}.                     T0​=Timu​.
初始速度                                              v                            0                                       \mathbf{v}_0                  v0​ 和偏置                                              b                            0                                       \mathbf{b}_0                  b0​ 设为零向量:
                                                    v                               0                                      =                            [                            0                            ,                            0                            ,                            0                                       ]                               T                                      ,                                                b                               0                                      =                                       [                                                                                                     b                                              a                                                                                                                                                  b                                              ω                                                                                                ]                                      =                            0.                                  \mathbf{v}_0 = [0, 0, 0]^T, \quad \mathbf{b}_0 = \begin{bmatrix} \mathbf{b}_a \\ \mathbf{b}_\omega \end{bmatrix} = \mathbf{0}.                     v0​=[0,0,0]T,b0​=[ba​bω​​]=0.
初始优化图: 利用先验因子初始化优化图:
  1.        // 初始位姿
  2.         prevPose_ = lidarPose.compose(lidar2Imu);  // 将激光雷达位姿转换为IMU位姿
  3.         gtsam::PriorFactor<gtsam::Pose3> priorPose(X(0), prevPose_, priorPoseNoise);  // 添加初始位姿先验因子
  4.         graphFactors.add(priorPose);
  5.         // 初始速度
  6.         prevVel_ = gtsam::Vector3(0, 0, 0);  // 初始速度为零
  7.         gtsam::PriorFactor<gtsam::Vector3> priorVel(V(0), prevVel_, priorVelNoise);  // 添加初始速度先验因子
  8.         graphFactors.add(priorVel);
  9.         // 初始偏差
  10.         prevBias_ = gtsam::imuBias::ConstantBias();  // 初始偏差为零
  11.         gtsam::PriorFactor<gtsam::imuBias::ConstantBias> priorBias(B(0), prevBias_, priorBiasNoise);  // 添加初始偏差先验因子
  12.         graphFactors.add(priorBias);
  13.         // 添加初始值
  14.         graphValues.insert(X(0), prevPose_);
  15.         graphValues.insert(V(0), prevVel_);
  16.         graphValues.insert(B(0), prevBias_);
复制代码

3. IMU积分与因子添加
IMU预积分,将IMU丈量值                                              a                            imu                                  ,                                   ω                            imu                                       \mathbf{a}_{\text{imu}}, \boldsymbol{\omega}_{\text{imu}}                  aimu​,ωimu​ 进行预积分,得到状态更新:
                                                                                           v                                                       k                                           +                                           1                                                                                                                                     =                                                       v                                           k                                                      +                                        (                                                       R                                           k                                                      ⋅                                        (                                                       a                                           imu                                                      −                                                       b                                           a                                                      )                                        +                                        g                                        )                                        ⋅                                        Δ                                        t                                        ,                                                                                                                            p                                                       k                                           +                                           1                                                                                                                                     =                                                       p                                           k                                                      +                                                       v                                           k                                                      ⋅                                        Δ                                        t                                        +                                                       1                                           2                                                      (                                                       R                                           k                                                      ⋅                                        (                                                       a                                           imu                                                      −                                                       b                                           a                                                      )                                        +                                        g                                        )                                        ⋅                                        Δ                                                       t                                           2                                                      ,                                                                                                                            R                                                       k                                           +                                           1                                                                                                                                     =                                                       R                                           k                                                      ⋅                                        exp                                        ⁡                                        (                                        (                                                       ω                                           imu                                                      −                                                       b                                           ω                                                      )                                        ⋅                                        Δ                                        t                                        )                                        .                                                                                \begin{aligned} \mathbf{v}_{k+1} &= \mathbf{v}_k + (\mathbf{R}_k \cdot (\mathbf{a}_{\text{imu}} - \mathbf{b}_a) + \mathbf{g}) \cdot \Delta t, \\ \mathbf{p}_{k+1} &= \mathbf{p}_k + \mathbf{v}_k \cdot \Delta t + \frac{1}{2} (\mathbf{R}_k \cdot (\mathbf{a}_{\text{imu}} - \mathbf{b}_a) + \mathbf{g}) \cdot \Delta t^2, \\ \mathbf{R}_{k+1} &= \mathbf{R}_k \cdot \exp((\boldsymbol{\omega}_{\text{imu}} - \mathbf{b}_\omega) \cdot \Delta t). \end{aligned}                     vk+1​pk+1​Rk+1​​=vk​+(Rk​⋅(aimu​−ba​)+g)⋅Δt,=pk​+vk​⋅Δt+21​(Rk​⋅(aimu​−ba​)+g)⋅Δt2,=Rk​⋅exp((ωimu​−bω​)⋅Δt).​
基于预积分丈量值构建:
                                         ImuFactor                            (                                       T                               k                                      ,                                       v                               k                                      ,                                       T                                           k                                  +                                  1                                                 ,                                       v                                           k                                  +                                  1                                                 ,                                       b                               k                                      ,                            PreintegratedIMU                            )                            .                                  \text{ImuFactor}(\mathbf{T}_{k}, \mathbf{v}_k, \mathbf{T}_{k+1}, \mathbf{v}_{k+1}, \mathbf{b}_k, \text{PreintegratedIMU}).                     ImuFactor(Tk​,vk​,Tk+1​,vk+1​,bk​,PreintegratedIMU).

  1.     const gtsam::PreintegratedImuMeasurements& preint_imu = dynamic_cast<const gtsam::PreintegratedImuMeasurements&>(*imuIntegratorOpt_);
  2.     gtsam::ImuFactor imu_factor(X(key - 1), V(key - 1), X(key), V(key), B(key - 1), preint_imu);
  3.     graphFactors.add(imu_factor);
  4.     // 添加IMU偏差之间的因子
  5.     graphFactors.add(gtsam::BetweenFactor<gtsam::imuBias::ConstantBias>(B(key - 1), B(key), gtsam::imuBias::ConstantBias(),
  6.                          gtsam::noiseModel::Diagonal::Sigmas(sqrt(imuIntegratorOpt_->deltaTij()) * noiseModelBetweenBias)));
复制代码
4. 里程计因子添加与优化
位姿因子,激光里程计位姿转换为IMU位姿后,构造先验因子:
                                         PriorFactor                            (                                       T                                           k                                  +                                  1                                                 ,                                       T                               cur                                      ,                                       Q                               pose                                      )                            ,                                  \text{PriorFactor}(\mathbf{T}_{k+1}, \mathbf{T}_{\text{cur}}, \mathbf{Q}_{\text{pose}}),                     PriorFactor(Tk+1​,Tcur​,Qpose​),
此中                                              T                            cur                                       \mathbf{T}_{\text{cur}}                  Tcur​ 为当前里程计位姿。
  1.     // 添加位姿因子
  2.     gtsam::Pose3 curPose = lidarPose.compose(lidar2Imu);
  3.     gtsam::PriorFactor<gtsam::Pose3> pose_factor(X(key), curPose, degenerate ? correctionNoise2 : correctionNoise);
  4.     graphFactors.add(pose_factor);
复制代码
优化求解,构建优化问题:
                                                    x                               ∗                                      =                            arg                            ⁡                                                   min                                  ⁡                                          x                                      ∑                            Factors                            .                                  \mathbf{x}^* = \arg\min_{\mathbf{x}} \sum \text{Factors}.                     x∗=argxmin​∑Factors.
调用优化器 optimizer.update 进行求解。
  1.     // 优化
  2.     optimizer.update(graphFactors, graphValues);
  3.     optimizer.update();
  4.     graphFactors.resize(0);
  5.     graphValues.clear();
复制代码

5. 优化后更新状态
优化完成后,从效果中提取最新状态:
                                                                                           T                                                       k                                           +                                           1                                                                                                                                     =                                        result                                        (                                                       X                                                           k                                              +                                              1                                                                     )                                        ,                                                                                                                            v                                                       k                                           +                                           1                                                                                                                                     =                                        result                                        (                                                       V                                                           k                                              +                                              1                                                                     )                                        ,                                                                                                                            b                                                       k                                           +                                           1                                                                                                                                     =                                        result                                        (                                                       B                                                           k                                              +                                              1                                                                     )                                        .                                                                                \begin{aligned} \mathbf{T}_{k+1} &= \text{result}(X_{k+1}), \\ \mathbf{v}_{k+1} &= \text{result}(V_{k+1}), \\ \mathbf{b}_{k+1} &= \text{result}(B_{k+1}). \end{aligned}                     Tk+1​vk+1​bk+1​​=result(Xk+1​),=result(Vk+1​),=result(Bk+1​).​
并重置IMU预积分器:
                                         imuIntegratorOpt                            −                            >                            resetIntegrationAndSetBias                            (                                       b                                           k                                  +                                  1                                                 )                            .                                  \text{imuIntegratorOpt}->\text{resetIntegrationAndSetBias}(\mathbf{b}_{k+1}).                     imuIntegratorOpt−>resetIntegrationAndSetBias(bk+1​).
  1.     // 为下一步重置预积分
  2.     gtsam::Values result = optimizer.calculateEstimate();
  3.     prevPose_  = result.at<gtsam::Pose3>(X(key));
  4.     prevVel_   = result.at<gtsam::Vector3>(V(key));
  5.     prevState_ = gtsam::NavState(prevPose_, prevVel_);
  6.     prevBias_  = result.at<gtsam::imuBias::ConstantBias>(B(key));
  7.     // 重置优化预积分对象
  8.     imuIntegratorOpt_->resetIntegrationAndSetBias(prevBias_);
  9.     // 检查优化是否失败
复制代码
6.重新传播
在重新传播之前,需要移除时间戳早于当前校正时间                                              t                            currentCorrectionTime                                  −                         Δ                         t                              t_{\text{currentCorrectionTime}} - \Delta t                  tcurrentCorrectionTime​−Δt 的IMU消息:
                                                    t                               imu                                      <                                       t                               currentCorrectionTime                                      −                            Δ                            t                            .                                  t_{\text{imu}} < t_{\text{currentCorrectionTime}} - \Delta t.                     timu​<tcurrentCorrectionTime​−Δt.
通过循环移除队列中满足上述条件的IMU消息,并记录最后一条被移除的IMU时间戳                                              t                            lastImuQT                                       t_{\text{lastImuQT}}                  tlastImuQT​:
                                                    t                               lastImuQT                                      =                                       t                               imu                                      .                                  t_{\text{lastImuQT}} = t_{\text{imu}}.                     tlastImuQT​=timu​.
  1.         // 首先弹出比当前校正数据旧的IMU消息
  2.         double lastImuQT = -1;
  3.         while (!imuQueImu.empty() && ROS_TIME(&imuQueImu.front()) < currentCorrectionTime - delta_t)
  4.         {
  5.             lastImuQT = ROS_TIME(&imuQueImu.front());
  6.             imuQueImu.pop_front();
  7.         }
复制代码
重新传播IMU
在IMU重新传播前,利用最新优化得到的偏置                                              b                            odom                                       \mathbf{b}_{\text{odom}}                  bodom​ 对IMU积分器重置:
  1. // 使用新优化的偏差重置偏差
  2. imuIntegratorImu_->resetIntegrationAndSetBias(prevBiasOdom);
复制代码
对于队列中剩余的IMU消息,从                                              t                            lastImuQT                                       t_{\text{lastImuQT}}                  tlastImuQT​ 开始,逐条进行积分:
                                                                                           v                                                       k                                           +                                           1                                                                                                                                     =                                                       v                                           k                                                      +                                        (                                                       R                                           k                                                      ⋅                                        (                                                       a                                           imu                                                      −                                                       b                                           a                                                      )                                        +                                        g                                        )                                        ⋅                                        Δ                                        t                                        ,                                                                                                                            p                                                       k                                           +                                           1                                                                                                                                     =                                                       p                                           k                                                      +                                                       v                                           k                                                      ⋅                                        Δ                                        t                                        +                                                       1                                           2                                                      (                                                       R                                           k                                                      ⋅                                        (                                                       a                                           imu                                                      −                                                       b                                           a                                                      )                                        +                                        g                                        )                                        ⋅                                        Δ                                                       t                                           2                                                      ,                                                                                                                            R                                                       k                                           +                                           1                                                                                                                                     =                                                       R                                           k                                                      ⋅                                        exp                                        ⁡                                        (                                        (                                                       ω                                           imu                                                      −                                                       b                                           ω                                                      )                                        ⋅                                        Δ                                        t                                        )                                        .                                                                                \begin{aligned} \mathbf{v}_{k+1} &= \mathbf{v}_k + (\mathbf{R}_k \cdot (\mathbf{a}_{\text{imu}} - \mathbf{b}_a) + \mathbf{g}) \cdot \Delta t, \\ \mathbf{p}_{k+1} &= \mathbf{p}_k + \mathbf{v}_k \cdot \Delta t + \frac{1}{2} (\mathbf{R}_k \cdot (\mathbf{a}_{\text{imu}} - \mathbf{b}_a) + \mathbf{g}) \cdot \Delta t^2, \\ \mathbf{R}_{k+1} &= \mathbf{R}_k \cdot \exp((\boldsymbol{\omega}_{\text{imu}} - \mathbf{b}_\omega) \cdot \Delta t). \end{aligned}                     vk+1​pk+1​Rk+1​​=vk​+(Rk​⋅(aimu​−ba​)+g)⋅Δt,=pk​+vk​⋅Δt+21​(Rk​⋅(aimu​−ba​)+g)⋅Δt2,=Rk​⋅exp((ωimu​−bω​)⋅Δt).​
此中:

  1.   // 从本次优化的开始处积分IMU消息
  2.             for (int i = 0; i < (int)imuQueImu.size(); ++i)
  3.             {
  4.                 sensor_msgs::Imu *thisImu = &imuQueImu[i];
  5.                 double imuTime = ROS_TIME(thisImu);
  6.                 double dt = (lastImuQT < 0) ? (1.0 / 500.0) : (imuTime - lastImuQT);
  7.                 imuIntegratorImu_->integrateMeasurement(gtsam::Vector3(thisImu->linear_acceleration.x, thisImu->linear_acceleration.y, thisImu->linear_acceleration.z),
  8.                                                         gtsam::Vector3(thisImu->angular_velocity.x,    thisImu->angular_velocity.y,    thisImu->angular_velocity.z), dt);
  9.                 lastImuQT = imuTime;
  10.             }
复制代码
免责声明:如果侵犯了您的权益,请联系站长,我们会及时删除侵权内容,谢谢合作!更多信息从访问主页:qidao123.com:ToB企服之家,中国第一个企服评测及商务社交产业平台。




欢迎光临 IT评测·应用市场-qidao123.com技术社区 (https://dis.qidao123.com/) Powered by Discuz! X3.4