VINS-mono代码阅读 -- 后端优化

tech2025-05-01  1

书接上回,初始化完成以后,就是进行非线性优化了,我们来看看这个函数吧。

1. 非线性优化

1) 状态向量

状态向量包括滑窗内的n+1个所有相机的状态,包括位置,朝向,速度,加速度计bias和陀螺仪bias,相机到IMU的外参,m+1个3D点的逆深度: X = [ x 0 , x 1 , ⋯   , x c b , λ 0 , λ 1 , ⋯   , λ m ] X = [x_0,x_1,\cdots,x_c^b,\lambda_0,\lambda_1,\cdots,\lambda_m] X=[x0,x1,,xcb,λ0,λ1,,λm] x k = [ p b k w , v b k w , q b k w , b a , b g ] x_k=[p_{b_k}^w,v_{b_k}^w,q_{b_k}^w,b_a,b_g] xk=[pbkw,vbkw,qbkw,ba,bg] x c b = [ p c b , q c b ] x_c^b=[p_c^b,q_c^b] xcb=[pcb,qcb]

2)目标函数

min ⁡ X = { ∥ r p − J p X ∥ 2 + ∑ k ∈ B ∥ r B ( z ^ b k + 1 b k , X ) ∥ P b k + 1 b k 2 + ∑ ( i , j ) ∈ C ∥ r C ( z ^ l C j , X ) ∥ p l C j 2 } \min_X=\lbrace \Vert r_p-JpX\Vert^2+\sum_{k\in B}\Vert r_B(\hat z_{b_{k+1}}^{b_k},X)\Vert _{P_{b_{k+1}}^{b_k}}^2 +\sum _{(i,j)\in C}\Vert r_C(\hat z_l^{C_j},X) \Vert_{p_l^{C_j}}^2\rbrace Xmin={rpJpX2+kBrB(z^bk+1bk,X)Pbk+1bk2+(i,j)CrC(z^lCj,X)plCj2} 三个残差项分别为:边缘化的先验信息,IMU的测量残差,视觉的重投影残差,三种残差都是用马氏距离表示的。

3)IMU约束

残差:两帧之间PVQ和bias的变化量 优化变量:

4)视觉约束

残差: 优化变量:

5)边缘化

根据次新帧是否是关键帧,分为两种边缘化策略:通过对比次新帧和次次新帧的视差,来决定marg掉次新帧或者最老帧,

当次新帧为关键帧时,MARGIN_OLD,将marg掉最老帧,及其看到的路标点和相关联的IMU数据,将其转换为先验信息加到整体的目标函数。当次新帧不是关键帧的时候,MARGIN_SECOND_NEW, 我们将直接丢掉次新帧及它的视觉观测边,而不对次新帧进行marg,因为我们认为当前帧和次新帧很相似,直接丢弃不会造成整个约束关系丢失过多信息,但是,要保留次新帧的IMU数据,从而保证IMU预积分的连贯性。

这里仅仅写出了优化的变量和对应的残差,具体的推导和代码的对应,需要单独的写一章出来。

2. 非线性优化的代码

1) 声明和引入鲁棒核函数

//1. 构建最小二乘问题 ceres::Problem problem; //2. 创建LossFunction对象, lossfunction 用来减少outlier的影响 ceres::LossFunction *loss_function; //loss_function = new ceres::HuberLoss(1.0); loss_function = new ceres::CauchyLoss(1.0);

2) 添加要优化的变量 位姿优化量

//3. 添加优化变量参数块 add vertex of 1.pose 2.speed 3.bias //添加要优化的变量:相机位姿、速度、加速度偏差、陀螺仪偏差 for (int i = 0; i < WINDOW_SIZE + 1; i++) { //定义了四元数的雅可比加法 ceres::LocalParameterization *local_parameterization = new PoseLocalParameterization(); //para_Pose[i]中存的是 第i帧的位姿,para_Pose[i]的大小SIZE_POSE是7 problem.AddParameterBlock(para_Pose[i], SIZE_POSE, local_parameterization); //SIZE_SPEEDBIAS 值为9 problem.AddParameterBlock(para_SpeedBias[i], SIZE_SPEEDBIAS); }

3)添加要优化的变量 相机到IMU的外参

如果不需要标定,就固定为常量

for (int i = 0; i < NUM_OF_CAM; i++) { ceres::LocalParameterization *local_parameterization = new PoseLocalParameterization(); problem.AddParameterBlock(para_Ex_Pose[i], SIZE_POSE, local_parameterization); if (!ESTIMATE_EXTRINSIC) { ROS_DEBUG("fix extinsic param"); problem.SetParameterBlockConstant(para_Ex_Pose[i]);//这个变量固定为常量 } else ROS_DEBUG("estimate extinsic param"); }

4)添加要优化的变量 Td

//添加要优化的变量 td if (ESTIMATE_TD) { problem.AddParameterBlock(para_Td[0], 1); //problem.SetParameterBlockConstant(para_Td[0]); }

5) vector2double()

给parameterBlock赋值。

void Estimator::vector2double() { for (int i = 0; i <= WINDOW_SIZE; i++) { para_Pose[i][0] = Ps[i].x(); para_Pose[i][1] = Ps[i].y(); para_Pose[i][2] = Ps[i].z(); Quaterniond q{Rs[i]}; para_Pose[i][3] = q.x(); para_Pose[i][4] = q.y(); para_Pose[i][5] = q.z(); para_Pose[i][6] = q.w(); para_SpeedBias[i][0] = Vs[i].x(); para_SpeedBias[i][1] = Vs[i].y(); para_SpeedBias[i][2] = Vs[i].z(); para_SpeedBias[i][3] = Bas[i].x(); para_SpeedBias[i][4] = Bas[i].y(); para_SpeedBias[i][5] = Bas[i].z(); para_SpeedBias[i][6] = Bgs[i].x(); para_SpeedBias[i][7] = Bgs[i].y(); para_SpeedBias[i][8] = Bgs[i].z(); } for (int i = 0; i < NUM_OF_CAM; i++) { para_Ex_Pose[i][0] = tic[i].x(); para_Ex_Pose[i][1] = tic[i].y(); para_Ex_Pose[i][2] = tic[i].z(); Quaterniond q{ric[i]}; para_Ex_Pose[i][3] = q.x(); para_Ex_Pose[i][4] = q.y(); para_Ex_Pose[i][5] = q.z(); para_Ex_Pose[i][6] = q.w(); } VectorXd dep = f_manager.getDepthVector(); for (int i = 0; i < f_manager.getFeatureCount(); i++) para_Feature[i][0] = dep(i); if (ESTIMATE_TD) para_Td[0][0] = td; }

6)先验信息

//4. 添加残差块 /*** * 每一项都是一个factor,这是ceres的使用方法,创建一个类继承于ceres::CostFunction类 * 重写Evaluate()函数定义residual的计算形式 */ //添加边缘化的残差信息,第一次进行优化的时候 last_marginalization_info 还为空值 if (last_marginalization_info) { // construct new marginlization_factor /*AddResidualBlock函数中第一个参数是 costfunction,第二个参数是lossfunction,第三个参数是参数快*/ MarginalizationFactor *marginalization_factor = new MarginalizationFactor(last_marginalization_info); problem.AddResidualBlock(marginalization_factor, NULL, last_marginalization_parameter_blocks); }

第一次执行这段代码的时候,没有先验信息,所以这段是不执行的,第二次执行的时候就有了,last_marginalization_info的赋值出现在后面的代码里。

class MarginalizationInfo { public: ~MarginalizationInfo(); int localSize(int size) const; int globalSize(int size) const; //加残差块相关信息(优化变量,待marg的变量) void addResidualBlockInfo(ResidualBlockInfo *residual_block_info); //计算每个残差对应的jacobian,并更新parameter_block_data void preMarginalize(); //pos为所有变量的维度,m要marg的变量,n为要保存的变量 void marginalize(); std::vector<double *> getParameterBlocks(std::unordered_map<long, double *> &addr_shift); //所有观测项 std::vector<ResidualBlockInfo *> factors; //m为要marg掉的变量个数,也就是parameter_block_idx的总localSize,以double为单位,VBias为9,PQ为6 //n为要保留下的优化变量的变量个数 n=localSize(parameter_block_size) - m int m, n; //<优化变量的内存地址, localSize> 为每个优化变量的变量大小,IMU(7,9,7,9) std::unordered_map<long, int> parameter_block_size; //global size int sum_block_size; //<待marg的优化变量内存地址,在parameter_block_index中的id> std::unordered_map<long, int> parameter_block_idx; //local size //<优化变量内存地址,数据> 每个优化块的数据 std::unordered_map<long, double *> parameter_block_data; std::vector<int> keep_block_size; //global size std::vector<int> keep_block_idx; //local size std::vector<double *> keep_block_data; //边缘化之后从信息矩阵中恢复出来的雅可比矩阵和残差向量 Eigen::MatrixXd linearized_jacobians; Eigen::VectorXd linearized_residuals; const double eps = 1e-8; };

当marg掉最老帧的时候,parameter_block_size为所有变量及其对应的localSize,parameter_block_data为对应的数据(double* 类型的)。下面举例说明,假设最老帧即第0帧观测到了k = 68个路标点那么代码中的parameter_block_size和parameter_block_data两个vector的大小均为: 2 ( V 0 : 1 ) + 11 ( P 0 : 10 ) + 1 ( T b c ) + 1 ( t b ) + 68 ( λ 1 : k ) = 83 2(V_{0:1})+11(P_{0:10})+1(T_{bc})+1(t_b)+68(\lambda_{1:k})=83 2(V0:1)+11(P0:10)+1(Tbc)+1(tb)+68(λ1:k)=83 pose=所有变量总localSize为: 2 × 9 + 11 × 6 + 1 × 6 + 1 × 1 + 68 × 1 = 159 2 \times 9+11\times6+1\times 6+1\times 1+68\times 1=159 2×9+11×6+1×6+1×1+68×1=159 parameter_block_idx为需要marg掉的变量,仍然假设第一帧看到了k=68个路标点,则parameter_block_idx数组的大小为: p a r a m e t e r _ b l o c k _ i d x . s i z e ( ) = 1 ( P 0 ) + 1 ( V 0 ) + 68 ( λ 1 : k ) = 70 parameter\_block\_idx.size()=1(P_0)+1(V_0)+68(\lambda _{1:k})=70 parameter_block_idx.size()=1(P0)+1(V0)+68(λ1:k)=70 MarginalizationInfo::m表示需要marg变量总的loaclSize, 即为 1 × 6 + 1 × 9 + 68 × 1 = 83 1\times 6+ 1\times 9+68\times 1=83 1×6+1×9+68×1=83 那么,MarginalizationInfo::n表示需要保留的变量的总localSIze,为 n = p o s − m = 159 − 83 = 76 n=pos-m=159-83=76 n=posm=15983=76

7) IMU残差

//添加IMU测量值残差 for (int i = 0; i < WINDOW_SIZE; i++) { int j = i + 1; if (pre_integrations[j]->sum_dt > 10.0) continue; IMUFactor* imu_factor = new IMUFactor(pre_integrations[j]); //AssResidualBlock 函数中的第一个参数是 costfunction 第二个参数是 lossfunction 第三个参数是 参数块 problem.AddResidualBlock(imu_factor, NULL, para_Pose[i], para_SpeedBias[i], para_Pose[j], para_SpeedBias[j]); }

8) 重投影残差

for (auto &it_per_id : f_manager.feature) { it_per_id.used_num = it_per_id.feature_per_frame.size(); if (!(it_per_id.used_num >= 2 && it_per_id.start_frame < WINDOW_SIZE - 2)) continue; ++feature_index; int imu_i = it_per_id.start_frame, imu_j = imu_i - 1; //第一个观测到该特征的帧 对应的特征点坐标 Vector3d pts_i = it_per_id.feature_per_frame[0].point; //遍历能观测到该特征的每个帧 for (auto &it_per_frame : it_per_id.feature_per_frame) { imu_j++; if (imu_i == imu_j) { continue; } Vector3d pts_j = it_per_frame.point; //ESTIMATED_TD : online estimate time offset between camera and imu if (ESTIMATE_TD) { ProjectionTdFactor *f_td = new ProjectionTdFactor(pts_i, pts_j, it_per_id.feature_per_frame[0].velocity, it_per_frame.velocity, it_per_id.feature_per_frame[0].cur_td, it_per_frame.cur_td, it_per_id.feature_per_frame[0].uv.y(), it_per_frame.uv.y()); problem.AddResidualBlock(f_td, loss_function, para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index], para_Td[0]); } else { ProjectionFactor *f = new ProjectionFactor(pts_i, pts_j); problem.AddResidualBlock(f, loss_function, para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index]); } f_m_cnt++; } }

9)重定位残差

//重定位残差 if(relocalization_info) { //printf("set relocalization factor! \n"); ceres::LocalParameterization *local_parameterization = new PoseLocalParameterization(); problem.AddParameterBlock(relo_Pose, SIZE_POSE, local_parameterization); int retrive_feature_index = 0; int feature_index = -1; //遍历特征 for (auto &it_per_id : f_manager.feature) { it_per_id.used_num = it_per_id.feature_per_frame.size(); if (!(it_per_id.used_num >= 2 && it_per_id.start_frame < WINDOW_SIZE - 2)) continue; ++feature_index; //获取观测到该特征的起始帧 int start = it_per_id.start_frame; if(start <= relo_frame_local_index) { while((int)match_points[retrive_feature_index].z() < it_per_id.feature_id) { retrive_feature_index++; } if((int)match_points[retrive_feature_index].z() == it_per_id.feature_id) { Vector3d pts_j = Vector3d(match_points[retrive_feature_index].x(), match_points[retrive_feature_index].y(), 1.0); Vector3d pts_i = it_per_id.feature_per_frame[0].point; ProjectionFactor *f = new ProjectionFactor(pts_i, pts_j); problem.AddResidualBlock(f, loss_function, para_Pose[start], relo_Pose, para_Ex_Pose[0], para_Feature[feature_index]); retrive_feature_index++; } } } }

3 求解非线性优化

//创建优化求解器 ceres::Solver::Options options; //设定优化器的solver_type 类型 options.linear_solver_type = ceres::DENSE_SCHUR; //options.num_threads = 2; //设定优化使用的方法,这里使用的是置信域类优化算法中的 dogleg options.trust_region_strategy_type = ceres::DOGLEG; //迭代次数 options.max_num_iterations = NUM_ITERATIONS; //options.use_explicit_schur_complement = true; //options.minimizer_progress_to_stdout = true; //options.use_nonmonotonic_steps = true; if (marginalization_flag == MARGIN_OLD) options.max_solver_time_in_seconds = SOLVER_TIME * 4.0 / 5.0; else options.max_solver_time_in_seconds = SOLVER_TIME; TicToc t_solver; ceres::Solver::Summary summary;//优化信息 //9. 开始求解 ceres::Solve(options, &problem, &summary); //cout << summary.BriefReport() << endl; ROS_DEBUG("Iterations : %d", static_cast<int>(summary.iterations.size())); ROS_DEBUG("solver costs: %f", t_solver.toc()); //求解完成后,将数组转换成向量 double2vector();

4. 边缘化

在这里只进行了边缘化,没有求解,求解的过程放在了下一轮的优化中。

1) marginalization_flag == MARGIN_OLD

把上一次先验项中的残差项(尺寸为n)传递给当前的先验项,并从中去除需要丢弃的状态量; MarginalizationInfo *marginalization_info = new MarginalizationInfo(); vector2double(); if (last_marginalization_info) { vector<int> drop_set; for (int i = 0; i < static_cast<int>(last_marginalization_parameter_blocks.size()); i++) { if (last_marginalization_parameter_blocks[i] == para_Pose[0] || last_marginalization_parameter_blocks[i] == para_SpeedBias[0]) drop_set.push_back(i); } // construct new marginlization_factor MarginalizationFactor *marginalization_factor = new MarginalizationFactor(last_marginalization_info); ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(marginalization_factor, NULL, last_marginalization_parameter_blocks, drop_set); marginalization_info->addResidualBlockInfo(residual_block_info); } 将滑窗内第0帧和第一帧之间的IMU预积分因子(pre_integrations[1])放到marginalization_info中 { if (pre_integrations[1]->sum_dt < 10.0) { IMUFactor* imu_factor = new IMUFactor(pre_integrations[1]); ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(imu_factor, NULL, vector<double *>{para_Pose[0], para_SpeedBias[0], para_Pose[1], para_SpeedBias[1]}, vector<int>{0, 1}); marginalization_info->addResidualBlockInfo(residual_block_info); } } 挑选出第一次观测帧为第0帧的路标点,将对应的多组视觉观测放到marginalization_Info中。 { int feature_index = -1; //遍历特征 for (auto &it_per_id : f_manager.feature) { it_per_id.used_num = it_per_id.feature_per_frame.size(); if (!(it_per_id.used_num >= 2 && it_per_id.start_frame < WINDOW_SIZE - 2)) continue; ++feature_index; int imu_i = it_per_id.start_frame, imu_j = imu_i - 1; if (imu_i != 0) continue; //当前特征所在的第一个观测帧中观测的点坐标 Vector3d pts_i = it_per_id.feature_per_frame[0].point; //遍历当前特征的观测帧 for (auto &it_per_frame : it_per_id.feature_per_frame) { imu_j++; if (imu_i == imu_j) continue; Vector3d pts_j = it_per_frame.point; if (ESTIMATE_TD) { ProjectionTdFactor *f_td = new ProjectionTdFactor(pts_i, pts_j, it_per_id.feature_per_frame[0].velocity, it_per_frame.velocity, it_per_id.feature_per_frame[0].cur_td, it_per_frame.cur_td, it_per_id.feature_per_frame[0].uv.y(), it_per_frame.uv.y()); ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(f_td, loss_function, vector<double *>{para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index], para_Td[0]}, vector<int>{0, 3}); marginalization_info->addResidualBlockInfo(residual_block_info); } else { ProjectionFactor *f = new ProjectionFactor(pts_i, pts_j); ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(f, loss_function, vector<double *>{para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index]}, vector<int>{0, 3}); marginalization_info->addResidualBlockInfo(residual_block_info); } } } } marginalization_info->preMarginalize(): 得到每次IMU和视觉观测(cost_function)对应的参数块(parameter_blocks),雅可比矩阵(jacobians),残差值(residuals)marginalization_info->marginalize(): 多线程计算整个先验的参数块,雅可比矩阵和残差值。 TicToc t_pre_margin; marginalization_info->preMarginalize(); ROS_DEBUG("pre marginalization %f ms", t_pre_margin.toc()); TicToc t_margin; marginalization_info->marginalize(); ROS_DEBUG("marginalization %f ms", t_margin.toc()); std::unordered_map<long, double *> addr_shift; for (int i = 1; i <= WINDOW_SIZE; i++) { addr_shift[reinterpret_cast<long>(para_Pose[i])] = para_Pose[i - 1]; addr_shift[reinterpret_cast<long>(para_SpeedBias[i])] = para_SpeedBias[i - 1]; } for (int i = 0; i < NUM_OF_CAM; i++) addr_shift[reinterpret_cast<long>(para_Ex_Pose[i])] = para_Ex_Pose[i]; if (ESTIMATE_TD) { addr_shift[reinterpret_cast<long>(para_Td[0])] = para_Td[0]; } vector<double *> parameter_blocks = marginalization_info->getParameterBlocks(addr_shift); if (last_marginalization_info) delete last_marginalization_info; //更新 last_marginalization_info last_marginalization_info = marginalization_info; last_marginalization_parameter_blocks = parameter_blocks; }//end of if(marginalization_flag == MARGIN_OLD)

2)marginalization_flag == MARGIN_SECOND_NEW

直接扔掉次新帧和它的视觉观测边,不对次新帧进行marg。但是保存了次新帧的IMU数据,从而保证预积分的连贯性。

else { /** * 当次新帧不是关键帧时,直接剪切掉次新帧和它的视觉边(该帧和路标点之间的关联),而不对次新帧进行marginalize处理 * 但是要保留次新帧的IMU数据,从而保证IMU预积分的连续性,这样才能积分计算出下一帧的测量值 */ if (last_marginalization_info && std::count(std::begin(last_marginalization_parameter_blocks), std::end(last_marginalization_parameter_blocks), para_Pose[WINDOW_SIZE - 1])) { MarginalizationInfo *marginalization_info = new MarginalizationInfo(); vector2double(); if (last_marginalization_info) { vector<int> drop_set; for (int i = 0; i < static_cast<int>(last_marginalization_parameter_blocks.size()); i++) { ROS_ASSERT(last_marginalization_parameter_blocks[i] != para_SpeedBias[WINDOW_SIZE - 1]); if (last_marginalization_parameter_blocks[i] == para_Pose[WINDOW_SIZE - 1]) drop_set.push_back(i); } // construct new marginlization_factor MarginalizationFactor *marginalization_factor = new MarginalizationFactor(last_marginalization_info); ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(marginalization_factor, NULL, last_marginalization_parameter_blocks, drop_set); marginalization_info->addResidualBlockInfo(residual_block_info); } TicToc t_pre_margin; ROS_DEBUG("begin marginalization"); //计算每次IMU和视觉观测(cost_function)对应参数块(parameters_block),雅可比矩阵(jacobians),残差值(residuals) marginalization_info->preMarginalize(); ROS_DEBUG("end pre marginalization, %f ms", t_pre_margin.toc()); TicToc t_margin; ROS_DEBUG("begin marginalization"); //多线程计算整个先验项的参数块,雅可比矩阵和残差值 marginalization_info->marginalize(); ROS_DEBUG("end marginalization, %f ms", t_margin.toc()); std::unordered_map<long, double *> addr_shift; for (int i = 0; i <= WINDOW_SIZE; i++) { if (i == WINDOW_SIZE - 1) continue; else if (i == WINDOW_SIZE) { addr_shift[reinterpret_cast<long>(para_Pose[i])] = para_Pose[i - 1]; addr_shift[reinterpret_cast<long>(para_SpeedBias[i])] = para_SpeedBias[i - 1]; } else { addr_shift[reinterpret_cast<long>(para_Pose[i])] = para_Pose[i]; addr_shift[reinterpret_cast<long>(para_SpeedBias[i])] = para_SpeedBias[i]; } } for (int i = 0; i < NUM_OF_CAM; i++) addr_shift[reinterpret_cast<long>(para_Ex_Pose[i])] = para_Ex_Pose[i]; if (ESTIMATE_TD) { addr_shift[reinterpret_cast<long>(para_Td[0])] = para_Td[0]; } vector<double *> parameter_blocks = marginalization_info->getParameterBlocks(addr_shift); if (last_marginalization_info) delete last_marginalization_info; last_marginalization_info = marginalization_info; last_marginalization_parameter_blocks = parameter_blocks; } }

优化部分到此结束。

最新回复(0)