在计算机视觉和机器人领域三维重建与运动估计是核心任务而光束法平差无疑是其背后最关键的优化引擎。无论是经典的SFM流程还是SLAM系统的后端优化BA的精度与鲁棒性直接决定了最终结果的质量。然而传统BA算法对初始值极为敏感一个糟糕的初始解往往会导致优化陷入局部极小值甚至直接发散。你是否也曾为如何给BA提供一个“足够好”的初始位姿和三维点而头疼近年来“免初始化”的BA方法逐渐进入研究视野它承诺摆脱对精确初始值的依赖直接从粗略甚至随机的猜测开始也能收敛到正确的解。这听起来很美好但实际效果如何不同方法之间孰优孰劣本文将以一篇前沿研究《Initialization-Free Bundle Adjustment Revisited》为引深入剖析免初始化BA的核心思想、主流方法特别是基于目标空间误差和变量投影的技术并通过结构化的对比为你呈现一份详实的“避坑指南”与工程实践参考。1. 光束法平差与初始化难题为什么这是个问题在深入“免初始化”方案之前我们必须先理解传统BA为何如此依赖初始化。1.1 光束法平差基础回顾光束法平差本质上是一个大规模非线性最小二乘问题。它的目标是同时优化相机参数姿态和内部参数和三维场景点的位置使得所有观测到的二维图像特征点与其对应的三维点投影到图像上的位置之间的重投影误差最小。用数学公式简要表示其代价函数 [ E \sum_{i1}^{m} \sum_{j1}^{n} \rho \left( | \mathbf{x}_{ij} - \pi(\mathbf{P}_i, \mathbf{X}_j) |^2 \right) ] 其中( \mathbf{x}_{ij} ) 是第 ( i ) 个相机观测到的第 ( j ) 个三维点 ( \mathbf{X}_j ) 的二维图像坐标。( \pi ) 是相机投影函数包含了相机的内参和外参姿态 ( \mathbf{P}_i ) 。( \rho ) 是鲁棒核函数如Huber损失用于抑制外点的影响。优化此问题通常使用高斯-牛顿法或列文伯格-马夸尔特法。这些方法都是局部迭代优化算法其共同特点是需要在优化起点即初始值附近代价函数近似为一个二次型才能保证快速、正确地收敛到全局最优解附近。1.2 初始化依赖的根源非线性与非凸性重投影误差模型 ( \pi(\mathbf{P}, \mathbf{X}) ) 是一个复杂的非线性函数特别是当相机模型包含镜头畸变时。这导致整个BA的代价函数是高度非凸的存在许多局部极小值。局部优化器的局限性LM/GN算法像“下山”的盲人只感知当前位置的坡度梯度和地形曲率海森矩阵近似。如果起点在一个小山谷的斜坡上它就会顺着坡滑到这个山谷的底部局部极小而无法知道远处是否存在更深的峡谷全局最优。初始值的角色一个好的初始值就是将这个“盲人”放在全局最优解所在的那个大峡谷的边缘让他能顺利滑到底部。一个差的初始值则可能把他放在另一个错误的山谷里或者放在一个陡峭的悬崖边导致优化步长计算失败而发散。在实际的SFM/SLAM流水线中我们通过对极几何、PnP、三角化等步骤来生成这些初始值。然而这些前置步骤本身也受噪声、外点和数值稳定性的影响产生的初始值可能并不“足够好”尤其是在大基线、低纹理或高速运动场景下。2. 免初始化BA的核心思想换个思路解决非凸问题既然问题出在非凸的代价函数和局部优化器上那么免初始化BA的思路无外乎两种1改造代价函数使其更“凸”2改造优化方法使其能跳出局部极小。2.1 目标空间误差从图像平面回到三维空间重投影误差是在2D图像平面上度量的误差。而目标空间误差是一种替代的误差度量方式。基本思想对于一对匹配的图像点 ( \mathbf{x}_1 ) 和 ( \mathbf{x}_2 )以及其对应的相机光心 ( \mathbf{C}_1 ) 和 ( \mathbf{C}_2 )可以反向投影形成两条射线。理想情况下这两条射线应在三维空间相交于一点 ( \mathbf{X} )。目标空间误差度量的就是这两条射线之间的“距离”例如最近点之间的距离或者点到射线的垂距。为什么可能免初始化OSE将误差定义在了3D空间其几何意义更直接。一些研究表明在某些参数化下基于OSE的代价函数其“盆地”可能比基于重投影误差的代价函数更宽即对初始值不那么敏感。优化过程可能更容易进入全局最优的吸引域。典型方法一些工作直接使用射线间距离作为代价。另一些方法则采用三角化然后重投影的交替优化策略将点坐标的优化从BA主问题中分离出来。2.2 变量投影降维打击与凸松弛变量投影是一种经典的优化技巧用于处理可分离的非线性最小二乘问题。BA问题恰好是“可分离”的给定相机参数三维点坐标可以通过线性最小二乘三角化最优解出反之亦然。基本思想将原问题 ( \min_{\mathbf{P}, \mathbf{X}} E(\mathbf{P}, \mathbf{X}) ) 转化为 ( \min_{\mathbf{P}} \tilde{E}(\mathbf{P}) )其中 ( \tilde{E}(\mathbf{P}) \min_{\mathbf{X}} E(\mathbf{P}, \mathbf{X}) )。即先将三维点坐标 ( \mathbf{X} ) 视为相机参数 ( \mathbf{P} ) 的函数通过最优三角化然后只对 ( \mathbf{P} ) 进行优化。这极大地降低了优化变量的维度。为什么可能免初始化降维变量数量大幅减少优化问题的复杂度降低搜索空间变小。隐式滤波在每一步对 ( \mathbf{X} ) 求最优解的过程等价于用当前最好的相机参数来重新三角化点这本身是一种对点坐标的“修正”可能有助于将优化路径引导至更好的区域。与凸松弛结合在只优化相机姿态的子问题上研究者更容易施加各种凸松弛技巧例如半定规划、拉格朗日对偶从而理论上保证找到全局最优解至少是更好的初始解。著名的Shor松弛在相对位姿估计中就有应用。2.3 其他辅助策略除了上述两种核心思想实践中常结合其他策略来提升鲁棒性鲁棒核函数的精心选择Huber、Cauchy、Tukey等核函数可以抑制外点影响防止少数错误匹配将优化带偏。自适应调整核函数阈值也是一个研究方向。多尺度或由粗到精优化先在低分辨率图像或简化模型如忽略畸变上进行优化得到一个粗略解再将其作为初始值用于全分辨率、完整模型的优化。随机初始化与重启最简单粗暴的方法。从多个随机初始点开始运行BA选择代价函数最小的结果。虽然计算量大但在离线处理中不失为一种实用方案。3. 环境准备与实验设置为了客观对比不同方法的性能我们需要一个可复现的实验环境。以下设置基于常见的研究实践。3.1 软件与库依赖操作系统Ubuntu 20.04 LTS 或 Windows 10/11 (WSL2推荐)。编程语言Python 3.8 或 C 14/17。Python更适合快速原型验证C用于高性能最终实现。核心计算库NumPy/SciPy用于矩阵运算和线性代数。Eigen(C)高性能线性代数库。Ceres Solver谷歌开源的通用非线性优化库内置BA求解器支持各种损失函数和参数化是BA实验的黄金标准。g2o另一个流行的图优化库在SLAM社区广泛应用。OpenCV用于图像处理、特征提取、相机模型和基础几何计算如三角化、PnP。数据集使用标准公开数据集进行评测是关键。Bundler/Photo Tourism Datasets如Notre Dame, Alamo, Trafalgar。ETH3D提供高精度地面真值适合定量评估。TUM RGB-D或EuRoC MAV针对SLAM的序列数据。自己生成仿真数据使用Blender或Unity等工具生成完全可控的3D场景和相机轨迹可以系统性地测试噪声、外点率和初始化偏差的影响。3.2 实验方案设计要点基准方法以Ceres或g2o实现的标准BA使用重投影误差LM优化作为基准。为其提供真实值添加噪声后的初始值模拟较好的初始化和随机扰动较大的初始值模拟差初始化。对比方法基于OSE的BA实现目标空间误差代价函数集成到Ceres中。变量投影法实现一个两阶段优化外层循环优化相机姿态使用SE(3)流形内层循环用线性最小二乘或SVD三角化所有点。鲁棒核函数调参对比不同核函数Huber vs Cauchy及阈值对收敛性的影响。评价指标收敛成功率在N次随机初始化的实验中成功收敛到“可接受解”的比例。如何定义“可接受”可以设定一个重投影误差阈值或者与地面真值的平均姿态/位置误差阈值。最终精度成功案例的最终重投影误差中位数/均值以及与地面真值的绝对轨迹误差。收敛速度迭代次数和计算时间。对初始偏差的鲁棒性系统性地增加初始位姿和点云的噪声水平观察各方法性能下降的曲线。4. 实战使用Ceres Solver对比标准BA与OSE BA下面我们通过一个简化的实例演示如何在Ceres Solver中实现并对比标准重投影误差和一种形式的目标空间误差。4.1 项目结构与依赖假设我们有一个简单的项目结构如下init_free_ba_study/ ├── data/ │ ├── cameras.txt # 相机初始参数带噪声 │ ├── points3d.txt # 3D点初始坐标带噪声 │ └── observations.txt # 2D观测数据 ├── include/ │ └── ba_cost_functions.h ├── src/ │ ├── standard_ba.cpp # 标准BA实现 │ ├── ose_ba.cpp # OSE BA实现 │ └── main.cpp # 主程序组织实验 ├── CMakeLists.txt └── README.mdCMakeLists.txt关键配置cmake_minimum_required(VERSION 3.10) project(InitFreeBA) set(CMAKE_CXX_STANDARD 14) find_package(Ceres REQUIRED) find_package(OpenCV REQUIRED) find_package(Eigen3 REQUIRED) include_directories(${Ceres_INCLUDE_DIRS} ${OpenCV_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIR}) add_executable(compare_ba src/main.cpp src/standard_ba.cpp src/ose_ba.cpp) target_link_libraries(compare_ba ${Ceres_LIBRARIES} ${OpenCV_LIBS})4.2 标准重投影误差代价函数实现include/ba_cost_functions.h:#ifndef BA_COST_FUNCTIONS_H #define BA_COST_FUNCTIONS_H #include ceres/ceres.h #include ceres/rotation.h #include Eigen/Core // 标准重投影误差代价函数 struct StandardReprojectionError { StandardReprojectionError(double observed_x, double observed_y, double fx, double fy, double cx, double cy) : observed_x(observed_x), observed_y(observed_y), fx(fx), fy(fy), cx(cx), cy(cy) {} template typename T bool operator()(const T* const camera_rotation, // 旋转向量 (角轴)3维 const T* const camera_translation, // 平移向量3维 const T* const point_3d, // 3D点坐标3维 T* residuals) const { // 1. 将点从世界坐标系转换到相机坐标系 T p[3]; ceres::AngleAxisRotatePoint(camera_rotation, point_3d, p); p[0] camera_translation[0]; p[1] camera_translation[1]; p[2] camera_translation[2]; // 2. 投影到归一化平面 T xp p[0] / p[2]; T yp p[1] / p[2]; // 3. 应用内参得到像素坐标 T predicted_x fx * xp cx; T predicted_y fy * yp cy; // 4. 计算残差 residuals[0] predicted_x - T(observed_x); residuals[1] predicted_y - T(observed_y); return true; } static ceres::CostFunction* Create(const double observed_x, const double observed_y, const double fx, const double fy, const double cx, const double cy) { return (new ceres::AutoDiffCostFunctionStandardReprojectionError, 2, 3, 3, 3( new StandardReprojectionError(observed_x, observed_y, fx, fy, cx, cy))); } private: double observed_x, observed_y; double fx, fy, cx, cy; }; #endif // BA_COST_FUNCTIONS_H4.3 目标空间误差代价函数实现以点到射线距离为例在include/ba_cost_functions.h中添加// 目标空间误差点到观测射线的距离。 // 注意此误差项仅依赖于相机姿态和3D点但需要另一个相机的观测来定义射线。 // 这里简化处理假设我们为每个3D点存储了其“参考观测”的相机光心方向。 struct ObjectSpaceError { ObjectSpaceError(const Eigen::Vector3d bearing_vector, const Eigen::Vector3d camera_center) : bearing_vector_(bearing_vector.normalized()), camera_center_(camera_center) {} template typename T bool operator()(const T* const point_3d, T* residuals) const { // bearing_vector_ 是单位化的观测方向在参考相机坐标系下但这里我们假设已转到世界系 // 更严谨的实现需要将bearing_vector_根据相机姿态进行旋转。 // 此处为演示做一个简化计算点到射线camera_center_ t * bearing_vector_的垂距。 Eigen::MatrixT, 3, 1 C camera_center_.castT(); Eigen::MatrixT, 3, 1 v bearing_vector_.castT(); Eigen::MatrixT, 3, 1 P(point_3d[0], point_3d[1], point_3d[2]); Eigen::MatrixT, 3, 1 P_C P - C; // 点到射线的距离 | (P-C) - ((P-C)·v) * v | T dot P_C.dot(v); Eigen::MatrixT, 3, 1 parallel dot * v; Eigen::MatrixT, 3, 1 perpendicular P_C - parallel; residuals[0] perpendicular.norm(); // 一个残差 return true; } // 注意这个CostFunction只优化3D点。在实际的BA中我们需要将相机姿态也作为参数 // 并建立姿态与bearing_vector的关系。这通常需要更复杂的参数化和多相机联合优化。 static ceres::CostFunction* Create(const Eigen::Vector3d bearing_vector, const Eigen::Vector3d camera_center) { return (new ceres::AutoDiffCostFunctionObjectSpaceError, 1, 3( new ObjectSpaceError(bearing_vector, camera_center))); } private: Eigen::Vector3d bearing_vector_; // 单位方向向量 Eigen::Vector3d camera_center_; };重要说明上述OSE实现是高度简化的。完整的OSE BA需要为每一对匹配的观测构建代价同时优化两个相机的姿态和一个3D点误差是两条射线之间的几何距离。实现复杂度较高此处仅为示意原理。4.4 主程序与实验逻辑src/main.cpp核心部分#include iostream #include vector #include random #include ba_cost_functions.h #include data_loader.h // 假设有一个加载数据的头文件 void RunStandardBA(const std::vectorObservation obs, std::vectorCamera cameras, std::vectorPoint3D points, double noise_level) { ceres::Problem problem; ceres::LossFunction* loss_function new ceres::HuberLoss(1.0); // 使用Huber核 for (const auto ob : obs) { Camera cam cameras[ob.camera_id]; Point3D pt points[ob.point_id]; ceres::CostFunction* cost_function StandardReprojectionError::Create(ob.x, ob.y, cam.fx, cam.fy, cam.cx, cam.cy); problem.AddResidualBlock(cost_function, loss_function, cam.rotation.data(), // 旋转向量指针 cam.translation.data(),// 平移向量指针 pt.data.data()); // 3D点指针 } // 设置参数化例如旋转使用角轴流形 for (auto cam : cameras) { problem.SetParameterization(cam.rotation.data(), new ceres::EigenQuaternionParameterization()); // 注意这里用四元数示例角轴有ceres::AngleAxisParameterization } ceres::Solver::Options options; options.linear_solver_type ceres::SPARSE_SCHUR; options.minimizer_progress_to_stdout true; options.max_num_iterations 100; ceres::Solver::Summary summary; ceres::Solve(options, problem, summary); std::cout Standard BA Report:\n summary.BriefReport() \n; } void AddNoiseToInitialGuess(std::vectorCamera cameras, std::vectorPoint3D points, double rot_noise_deg, double trans_noise, double point_noise) { std::default_random_engine generator; std::normal_distributiondouble rot_dist(0.0, rot_noise_deg * M_PI / 180.0); std::normal_distributiondouble trans_dist(0.0, trans_noise); std::normal_distributiondouble point_dist(0.0, point_noise); for (auto cam : cameras) { // 为旋转添加噪声需转换为合适的表示形式添加此处简化 // 为平移添加噪声 for (int i 0; i 3; i) { cam.translation[i] trans_dist(generator); } } for (auto pt : points) { for (int i 0; i 3; i) { pt.data[i] point_dist(generator); } } } int main() { // 1. 加载数据地面真值 auto [true_cameras, true_points, observations] LoadGroundTruthData(data/ground_truth.txt); // 2. 生成带噪声的初始值 std::vectorCamera init_cameras true_cameras; std::vectorPoint3D init_points true_points; double noise_level 10.0; // 噪声水平旋转10度平移和点坐标10%的尺度 AddNoiseToInitialGuess(init_cameras, init_points, noise_level, noise_level * 0.1, noise_level * 0.1); // 3. 运行标准BA std::cout Running Standard BA with Noisy Initialization \n; auto cameras_ba1 init_cameras; auto points_ba1 init_points; RunStandardBA(observations, cameras_ba1, points_ba1, noise_level); double error1 ComputeReprojectionError(cameras_ba1, points_ba1, observations); // 4. 运行OSE BA (简化版可能需要不同的数据结构和优化循环) // std::cout Running OSE BA (Simplified) \n; // RunOSEBA(...); // 5. 评估结果 std::cout Final Reprojection Error (Standard BA): error1 pixels\n; // 比较与地面真值的轨迹误差等... return 0; }4.5 运行与结果分析编译并运行程序后观察输出。在噪声水平较低时标准BA和OSE BA都应能成功收敛。逐步增大noise_level你会观察到标准BA在初始噪声超过一定阈值后优化可能开始发散最终误差极大或收敛到一个错误的局部极小值误差比地面真值对应的误差大很多。OSE BA其收敛曲线可能更加平缓在更大的初始噪声范围内仍能收敛到可接受的解。但最终精度可能略低于在良好初始化下的标准BA因为OSE误差本身是重投影误差的一种近似。关键输出指标summary.termination_type: 查看是收敛 (CONVERGENCE) 还是失败 (NO_CONVERGENCE,FAILURE)。summary.final_cost: 最终代价函数值。summary.iterations: 迭代次数。自定义计算的平均重投影误差和绝对轨迹误差。5. 常见问题与排查思路在实现和实验免初始化BA时你可能会遇到以下典型问题问题现象可能原因排查与解决思路优化不收敛代价函数值NaN或inf1. 初始值太差导致投影计算出现非法值如深度为负。2. 参数化错误特别是旋转表示四元数需单位化角轴需小量。3. 损失函数参数过于激进或观测数据中存在严重外点。1. 检查初始深度值确保点在相机前方。可尝试从随机初始值多次重启。2. 使用Ceres内置的流形如EigenQuaternionParameterization,AngleAxisParameterization。3. 调整鲁棒核函数的阈值或先进行严格的外点滤除如RANSAC。OSE BA收敛但最终重投影误差比标准BA大OSE误差是几何距离与像素级重投影误差最小化目标不完全等价。优化OSE得到的解在重投影误差度量下未必最优。这是预期之内的权衡。可以混合使用先用OSE BA得到一个较好的初始解再将其作为标准BA的输入进行精优化。变量投影法实现复杂且速度慢变量投影需要在外层循环优化姿态内层循环对所有点进行三角化。三角化步骤如果实现不当如直接求伪逆会成为性能瓶颈。1. 使用高效的线性三角化方法如SVD。2. 注意三角化也需要处理多视图情况可能涉及迭代重加权。3. 考虑使用舒尔补技巧Ceres的SPARSE_SCHUR求解器本质上利用了BA问题的稀疏结构与变量投影有思想关联。对极端外点敏感任何BA方法都怕外点。免初始化方法可能在糟糕的初始阶段更容易被外点带偏。1.前端是关键提升特征匹配质量使用更鲁棒的描述子和匹配器。2.RANSAC预处理在BA之前使用RANSAC进行几何验证对极几何、PnP。3.动态核函数使用自适应阈值的鲁棒核或在优化过程中逐步收紧阈值。6. 最佳实践与工程建议基于现有研究和实验以下是在实际项目中应用或选择BA初始化策略的建议没有银弹分层处理是王道离线SFM计算资源允许采用由粗到精的管道。先使用低分辨率图像、简化相机模型和OSE或变量投影法进行全局运动恢复结构得到一个粗糙但全局一致的解。再以此作为初始值进行完整的、高精度的标准BA优化。在线SLAM对实时性要求高。通常依赖短期跟踪的累积来提供较好的初始值如通过IMU预积分、视觉里程计。在这种情况下标准BA或其增量式变种如iSAM2已足够鲁棒。免初始化BA可能因计算量大而不适用。理解问题尺度与场景小尺度、相机运动平缓标准BA对初始化不敏感简单初始化如恒等矩阵或线性方法八点法足以。大尺度、存在闭环、初始估计可能很差这是免初始化BA最能发挥价值的场景。考虑将OSE BA或基于凸松弛的全局方法作为回环检测后的全局优化步骤。充分利用现有库不要重复造轮子Ceres Solver和g2o已经极度优化。你的主要工作应是定义正确的代价函数和参数化以及设计好的优化流程如多阶段优化而非实现优化算法本身。在Ceres中可以方便地切换损失函数、线性求解器并添加各种参数块使用LocalParameterization或Manifold。初始化的质量评估在运行完整BA之前实现一个快速检查。例如计算初始解下的平均重投影误差。如果误差巨大如数百像素那么标准BA大概率会失败应触发备用方案如降级到OSE BA或全局搜索。变量投影的工程实现技巧如果决定实现变量投影注意姿态参数化。使用李代数如Sophus库可以方便地在流形上进行优化。内层的点三角化可以并行化因为每个点的三角化是独立的。考虑使用子空间划分先优化一部分关键相机姿态和点固定其余的交替进行以降低单次优化问题的规模。仿真与定量评估不可或缺在将新方法应用到真实数据前必须在仿真数据上进行系统测试。仿真可以控制噪声水平、外点率、场景类型和初始偏差从而清晰地揭示方法的优缺点和失效边界。免初始化BA是一个活跃的研究领域它并非要完全取代传统BA而是为其在困难场景下提供更鲁棒的启动方式。对于工程师而言核心在于理解不同误差度量重投影误差 vs 目标空间误差和优化策略联合优化 vs 变量投影的利弊并根据具体应用场景数据规模、实时性要求、初始化质量构建一个混合、分层的优化管道。将OSE或全局方法作为“粗调”模块将标准BA作为“精调”模块往往是兼顾鲁棒性与精度的最佳实践。通过本文的代码示例和实验框架你可以开始自己的探索在三维视觉项目中更好地驾驭光束法平差这一强大工具。