增量式SFM自由网平差:从数学原理到工程实践
1. 从“拍照片”到“建模型”为什么我们需要增量式SFM自由网平差如果你玩过三维重建或者对计算机视觉、摄影测量有点兴趣那你肯定听过“SFM”Structure from Motion从运动恢复结构这个词。简单说就是给你一堆从不同角度拍摄的同一个物体的照片SFM能帮你算出每张照片是在什么位置、什么角度拍的这叫相机位姿以及这个物体上各个点的三维坐标这叫场景结构。听起来很酷对吧但当你真正上手去做尤其是处理几十、上百甚至上千张无序照片时问题就来了计算量爆炸、误差累积、结果漂移……最后重建出来的模型可能歪七扭八甚至直接散架。这时候“增量式SFM”和“自由网平差”Bundle Adjustment, BA就成了解决问题的核心搭档。增量式SFM就像搭积木从两张照片开始一点点把新的照片和三维点加进来逐步扩大重建规模。而BA特别是我们今天要聊的“无先验约束下的自由网平差”就是每次搭完一部分积木后进行一次全局的“微调”和“压实”确保整个结构又稳又准。这个“无先验约束”是关键意味着我们不依赖任何外部信息比如GPS坐标、已知的控制点来固定坐标系完全让数据自己“说话”通过数学优化找到一个内部最一致的三维模型。为什么这很重要因为现实世界中我们很多时候就是拿个手机或普通相机随手一拍根本没有精确的测量设备。自由网平差能让这些“野生”数据自己对齐生成一个比例正确、形状准确但可能整体位置和朝向是任意的的三维模型。这对于考古数字化、室内场景重建、甚至影视特效中的场景扫描都是基础且关键的一步。理解了它你才算真正摸到了SFM的“里子”而不是停留在调用OpenMVG或COLMAP的API层面。2. 核心基石最小二乘问题在SFM中的本质要搞懂自由网平差必须先吃透它的数学心脏——非线性最小二乘Non-linear Least Squares。在SFM的BA问题里我们到底在优化什么一句话让所有照片上观察到的二维像素点与我们根据当前估计的三维点和相机参数“投影”回去的二维点之间的差距最小。2.1 重投影误差一切优化的目标函数假设我们有一个三维点 (X_j [x, y, z]^T)在世界坐标系下和一个相机其位姿用旋转矩阵 (R_i) 和平移向量 (t_i) 表示内参矩阵为 (K_i)。这个三维点投影到第i张照片上的像素坐标应该是[ \hat{u}_{ij} \pi(K_i (R_i X_j t_i)) ]这里 (\pi) 是投影函数把相机坐标系下的三维点 ([X_c, Y_c, Z_c]^T) 变成归一化平面坐标 ([X_c/Z_c, Y_c/Z_c]^T)再乘以内参得到像素坐标 ([u, v]^T)。而在第i张照片上我们通过特征匹配实际观测到这个点位于 (u_{ij} [u_{obs}, v_{obs}]^T)。那么重投影误差就是[ e_{ij} u_{ij} - \hat{u}_{ij} ]BA的目标就是调整所有相机参数所有 (R_i, t_i, K_i)和所有三维点坐标所有 (X_j)使得所有这样的误差的平方和最小[ \min_{{R_i, t_i, K_i}, {X_j}} \sum_{i, j} \rho(| e_{ij} |^2) ]这里的 (\rho) 是一个鲁棒核函数比如Huber核用来降低误匹配外点的影响这不是我们今天讨论的重点但它是工程实践中不可或缺的一环。2.2 问题的“病态”与自由网的特殊性如果我们直接对上述目标进行优化会立刻遇到一个问题尺度模糊性。想象一下你把所有相机和三维点同时放大两倍投影到图像上的像素坐标是不变的。同样把整个场景旋转和平移一下只要所有相机和点一起动投影结果也不变。这意味着我们的优化问题有多个等价的解数学上称之为“零空间”Null Space或“自由度”Gauge Freedom。具体来说一个三维欧几里得空间中的SFM问题在没有先验信息时存在7个自由度3个平移、3个旋转、1个尺度。这导致问题的Hessian矩阵二阶导数矩阵是奇异的直接求逆会失败。这就是“自由网”需要特别处理的原因。我们需要以一种优雅的方式“固定”这个网消除这7个自由度让优化问题变得良态Well-posed同时又不能引入虚假的约束扭曲模型本身。3. 增量式SFM流程中的BA集成策略增量式SFM不是一上来就处理所有照片它有一个清晰的流水线。BA作为其中的优化模块被巧妙地集成在关键节点上。3.1 初始化构建第一个稳定的“种子模型”一切从两张图片开始。通常选择匹配点数量多、视差足够大保证三角化精度的一对图像。相对姿态估计通过五点法或八点法RANSAC计算这两张图片之间的本质矩阵E或基础矩阵F分解得到相对旋转 (R_{rel}) 和相对平移 (t_{rel})带尺度模糊。三角化利用这对相对姿态将匹配的特征点三角化成最初的一批三维点 (X_j)。此时我们通常将第一个相机的坐标系设为世界坐标系(R_1 I, t_1 0)第二个相机位姿即为估计出的相对姿态。第一次双视图BA立即对这两个相机和第一批三维点进行一次小规模的BA。这次BA已经是自由网平差但因为它只涉及两个相机固定第一个相机锁定其旋转平移为零和其中一个三维点的深度或固定平移向量的范数通常就能消除尺度模糊为整个重建确立初始尺度和坐标系。这是一个非常关键的步骤初始化的精度直接影响后续增量添加的稳定性。3.2 增量添加与部分BA初始化后进入循环下一视图选择从剩余图像中选择能看到当前已重建三维点数量最多的图片。这保证了新视图有足够的2D-3D对应关系用于姿态估计。姿态估计PnP利用当前三维点和新视图中的2D特征点通过EPnP或UPnP等算法配合RANSAC抗噪求解新相机的位姿 (R_{new}, t_{new})。三角化新点新视图加入后它与已有视图之间会产生新的匹配点对将这些点三角化增加三维点云的数量和密度。局部BAPartial BA这是增量式SFM保持效率的核心。我们不会每次都优化所有参数。通常只优化新加入的相机、与它相关的三维点即能被它看到的点以及最近加入的几个其他相机。其他远离当前添加区域的“旧”相机和点则保持固定。这极大地减少了每次BA的变量规模加快了计算速度。此时的平差仍然是在当前定义的局部坐标系下进行的自由网平差。3.3 全局BA与“浮动”的自由网当增量添加进行到一定程度例如每添加50张图后或者所有图像都添加完毕后必须执行一次全局BA。这次所有相机参数和所有三维点都参与优化。这才是“无先验约束下的自由网平差”发挥全部威力的舞台。此时我们面对的是成千上万个待优化参数。直接求解会触发前面提到的尺度、旋转、平移自由度问题。解决方法不是强行固定某个相机或点那会扭曲优化结果而是采用更数学化的方式使用先验权重或添加软约束一种常见实践是在目标函数中为选定的“基准”相机或点添加一个非常弱的先验项例如 (\lambda | t_1 |^2)其中 (\lambda) 是一个极小的数如 (10^{-8})。这相当于告诉优化器“理论上 (t_1) 应该是零但我并不强求只是稍微倾向于零”。这种方法能稳定数值解又不会对结果产生可观测的影响。另一种等价的实现方式是在求解线性方程如高斯-牛顿法中的正规方程时对对应的参数块在Hessian矩阵的对角线上加一个小的阻尼因子。经过全局BA优化后我们得到了一个在内部度量意义下最优的三维模型。它的形状、相对距离、角度都是准确的但它整体漂浮在一个任意的坐标系中可能旋转了、平移了、缩放了。对于很多应用如纹理映射、三维浏览这已经完全足够。4. 自由网平差的数值实现与关键技巧理论很美但落地到代码里全是细节。这里分享几个在实现自由网平差时直接影响结果好坏和速度快慢的关键点。4.1 参数化与雅可比矩阵计算优化需要计算重投影误差对每个参数的导数雅可比矩阵。参数化方式影响收敛速度和稳定性。三维点直接用欧几里得坐标 ([x, y, z]^T) 即可雅可比矩阵推导相对直接。相机旋转千万不要直接用3x3矩阵的9个元素这引入了6个冗余约束。主流做法是使用李代数 (so(3)) 上的三维向量 (\omega)旋转向量或四元数。在优化迭代中我们通常在小扰动模型下工作(R \leftarrow R \cdot \text{Exp}(\delta \omega))其中 (\text{Exp}) 是指数映射将旋转向量转为旋转矩阵。这样优化变量就是三维的 (\delta \omega)雅可比矩阵也是关于这个三维扰动的导数。相机平移直接用三维向量 (t)。相机内参对于针孔模型通常优化焦距 (f_x, f_y) 和主点 (c_x, c_y)。如果考虑畸变则加上径向畸变系数 (k_1, k_2, k_3) 等。雅可比矩阵的计算涉及链式法则需要耐心推导或利用自动微分库如Ceres Solver中的Jet类型。手工推导时一个常见的技巧是利用归一化平面坐标 (p [X_c/Z_c, Y_c/Z_c, 1]^T) 作为中间变量能简化计算。4.2 稀疏性与Schur Complement技巧BA问题的Hessian矩阵 (J^T J)高斯-牛顿法或近似Hessian列文伯格-马夸尔特法有一个美妙的稀疏块结构。误差项只与看到它的相机和它本身这个三维点有关。因此Hessian矩阵可以按相机块C和点块P进行分块[ H \begin{bmatrix} B E \ E^T C \end{bmatrix} ] 其中 (B) 是对角块矩阵相机-相机二阶导(C) 也是对角块矩阵点-点二阶导(E) 是连接相机和点的非对角块矩阵。在求解增量方程 (H \delta x -g) 时直接对巨大的 (H) 求逆不可行。利用Schur Complement技巧我们可以先消去三维点参数 (\delta X)通常点数远多于相机数得到一个只关于相机参数增量 (\delta C) 的、规模小得多的方程[ (B - E C^{-1} E^T) \delta C -g_C E C^{-1} g_P ]这个新的系数矩阵 (S B - E C^{-1} E^T) 被称为舒尔补矩阵它只与相机有关且保持了稀疏性相机i和j相连仅当它们看到同一个三维点。求解出 (\delta C) 后再回代求解 (\delta X)。这是所有高效BA库如Ceres, g2o的核心。在自由网平差中这个技巧同样适用我们只需要在处理那个奇异的、具有零空间的全局Hessian矩阵时格外小心。4.3 处理自由度的工程实践阻尼与基准选择在实际代码中如何处理那7个自由度我个人的经验是隐式处理推荐使用像Ceres Solver这样的成熟库。它内部在求解线性系统时会检测并处理数值奇异性。当你使用DENSE_SCHUR或SPARSE_SCHUR求解器时它通常能稳健地处理自由网问题。其原理类似于在奇异方向添加微小的阻尼使求解器能找到一个最小范数解。显式固定如果你自己实现求解器一个简单有效的方法是在全局BA迭代开始前选定一个相机比如第一个相机和一个三维点比如第一个被三角化的点。固定该相机的旋转为单位矩阵平移为零。固定该三维点的深度为一个常数或者固定其XY坐标为零Z为一个正数。在后续的每次线性方程求解中将这些固定参数对应的行和列从Hessian矩阵和梯度向量中直接移除。这相当于显式地消除了7个自由度。这种方法直观但要注意基准选择的任意性不会影响内部几何。注意固定基准时务必确保你固定的是“一个相机”和“一个点”的足够多的参数来消除所有自由度。一个常见的错误是只固定了相机的平移但旋转和尺度自由度依然存在。更稳妥的方法是固定第一个相机的全部6个参数旋转和平移这直接定义了世界坐标系的原点和朝向同时通过整个优化过程隐式地确定了尺度因为重投影误差是尺度不变的但优化过程会收敛到一个使误差最小的特定尺度上这个尺度由初始化决定。4.4 鲁棒核函数不让一个误匹配毁掉整个模型这是BA能否实用的关键。重投影误差服从高斯分布只是理想假设特征匹配中必然存在外点误匹配。如果直接用平方误差L2范数这些外点会产生巨大的误差把优化拉偏。必须使用鲁棒核函数例如Huber损失或Cauchy损失。它们的作用是当误差小于某个阈值时行为类似平方误差保证精度当误差大于阈值时增长变缓如变成线性增长从而抑制外点的影响。在Ceres中这只需要添加一个LossFunctionceres::Problem problem; ceres::LossFunction* loss_function new ceres::HuberLoss(1.0); // 阈值通常取像素单位如1.0 problem.AddResidualBlock(cost_function, loss_function, parameters);阈值的选择有讲究太小会过度抑制正常数据太大会放过外点。通常根据特征点定位精度来设比如0.5到2个像素。可以先用一个较大阈值如4.0进行几轮优化剔除误差大的点再用较小阈值如1.0进行精优化。5. 实战中的坑与性能优化心法理论跑通只是第一步让它在真实数据上稳定、高效地工作才是挑战。下面是我踩过的一些坑和总结的优化技巧。5.1 初始化失败与视图选择策略增量式SFM最脆弱的环节是初始化。如果初始的两视图重建质量差误差会像雪球一样越滚越大。坑1视差不足选择的初始图像对视差太小导致三角化的三维点深度不确定性极大Z值接近零甚至为负。对策在选择初始对时不仅要看匹配数量还要计算本质矩阵分解后三角化点的正向深度比例和视差角。确保大部分点有合理的深度和足够的视差通常大于5度。坑2纯旋转或近似纯旋转如果初始图像对之间主要是旋转没有足够的平移则无法三角化尺度也无法确立。对策在匹配点对中检查归一化坐标如果对极几何约束主要来自基础矩阵F而非本质矩阵E则可能是纯旋转场景需要避免将其作为初始对。视图选择策略也至关重要。贪心地选择“看到最多已有点”的视图并不总是最优因为这可能导致重建沿着一个“主干”延伸而忽略了其他视角使得后续BA难以纠正早期累积的漂移。一种改进策略是引入“探索性”偶尔选择能看到已重建区域边缘或未覆盖区域的图像增加模型的闭合环Loop这对BA纠正漂移有巨大帮助。5.2 尺度漂移与闭环检测在长序列或无GPS数据的重建中尺度漂移是自由网平差的顽疾。由于误差累积模型末端的尺度可能与起始端不一致。现象重建出的长廊尽头比实际更窄或更宽一个圆形轨迹重建出来变成了螺旋。根本原因增量过程中每次PnP和新点三角化都引入微小误差BA虽然能局部优化但在没有全局约束的长序列中无法完全消除跨远距离的累积误差。终极解药闭环检测Loop Closure。当系统识别出当前图像与很久以前添加的图像匹配成功时就形成了一个闭环。闭环提供了跨越长距离的强约束。在BA中闭环对应的重投影误差项将直接连接时间上相隔很远的相机从而给优化器提供了纠正漂移的“锚点”。实现闭环后执行全局BA通常能显著改善尺度一致性和整体精度。现代SFM系统如COLMAP的全局BA模块都会紧密集成闭环检测信息。5.3 数值稳定性与参数缩放BA优化中参数的量级可能差异巨大。旋转的扰动 (\delta \omega) 单位是弧度平移 (t) 单位可能是米三维点坐标 (X) 单位也是米内参焦距 (f) 单位是像素。如果直接把这些参数堆成一个向量Hessian矩阵的条件数会很大导致线性求解器数值不稳定收敛缓慢甚至失败。技巧参数缩放Parameter Scaling。在构造优化问题前对参数进行归一化。一个常用的做法是将整个场景的坐标进行归一化使得所有三维点的重心位于原点且点云的平均距离原点的RMS距离为1或一个固定值。这相当于对平移和点坐标进行了缩放。同时确保旋转扰动和焦距等参数也在合理的量级内。在Ceres中可以使用Problem::SetParameterBlockConstant和Manifold来隐式处理但更根本的是在数据预处理阶段做好尺度归一化。5.4 内存与速度优化从稠密到稀疏的抉择全局BA可能涉及上万个相机和百万个点雅可比矩阵和Hessian矩阵非常庞大。使用稀疏求解器如前所述利用舒尔补和Hessian的稀疏性是必须的。Ceres中的SPARSE_SCHUR求解器会使用SuiteSparse或Eigen的稀疏Cholesky分解效率远高于稠密求解。限制优化规模对于超大场景即使稀疏求解也吃力。可以采用“分层BA”或“子图BA”。将整个场景划分为多个子图先在子图内进行BA然后将子图视为一个“超级节点”再进行子图间的BA。或者只优化最近的一部分参数类似局部BA定期执行全局BA但使用更低的频率或更粗糙的采样。GPU加速BA中的大量运算是矩阵-向量乘和雅可比矩阵计算非常适合并行。一些库如Ceres支持使用CUDA进行部分计算的加速。对于实时SLAMGPU加速的BA几乎是标配。6. 进阶思考从自由网到有控网掌握了无先验约束的自由网平差你就拥有了从无序图像中“无中生有”构建三维世界的能力。但它的输出是一个“浮动”的模型。如何将它锚定到真实世界尺度恢复如果你知道场景中某两点的真实距离比如地砖是30厘米可以在BA后通过相似变换旋转、平移、缩放将模型对齐到这个已知距离上。绝对定向如果你有少量控制点GCPs的真实世界坐标如通过RTK测量可以在BA的目标函数中加入这些控制点的三维坐标误差项进行有控网平差。这会将模型牢固地绑定到特定坐标系如UTM坐标系下。融合IMU/GPS在移动端或无人机场景可以将BA与IMU惯性测量单元和GPS数据进行紧耦合或松耦合优化即视觉惯性里程计VIO或SLAM直接获得尺度、重力方向甚至绝对位置。自由网平差是这一切的基础。它确保了视觉观测的内部一致性为后续的融合与对齐提供了一个最优的、自洽的初始解。理解它的每一个细节从数学原理到工程实现从增量策略到调参技巧是打通SFM任督二脉的不二法门。下次当你看到COLMAP的进度条在“Running bundle adjustment...”时希望你能会心一笑知道那背后正是一场精妙绝伦的、关于几何与优化的交响乐。