简介：面向摄影测量与计算机视觉学习者的光束法区域网平差VC++实现，用于解算多张像片的外方位元素并优化特征点三维坐标，适合正在学习结构从运动、立体匹配或三维重建的开发者参考。资源包共32个文件，以8个h头文件、7个cpp源文件为主，辅以5个txt数据/结果文档及VC++工程配置（dsp/dsw/rc等），整体仅55KB，代码结构清晰，便于直接阅读与编译调试。已有1320人学习下载，说明其受到一定关注。内容涵盖光束法平差核心流程：基于最小二乘原理建立误差方程，利用Levenberg-Marquardt迭代优化，同时考虑像片旋转角与平移量；工程中通过矩阵类封装实现数值计算，并提供原始数据与结果文本，可对照验证外方位元素解算精度。对于想理解摄影测量平差原理、或需要迁移到无人机测绘/遥感应用的开发者，这份紧凑的代码包具有直接参考价值。
1. 光束法区域网平差：为什么它是摄影测量的最终裁判
很多人在学摄影测量时都有过这个困惑：单模型空间后方交会就能求出外方位元素，逐模型拼接也能出成果，为什么还要搞一套把所有影像捆在一起的光束法区域网平差？反直觉的地方在于，光束法把所有影像的所有光线全部放进同一个大型方程系统，同时解算每张影像的6个外方位元素和每个待定点的3个物方坐标。它不是局部打补丁，而是让成千上万条光束在最小二乘意义下同时穿过物方点，因此它在数学上是最严格、整体精度最高的一类解法。这份思路落到 VC++ 代码里，就是从一条共线方程到全套法方程求解的完整骨架。弄懂光束法平差的底层机制，才算真正打开了摄影测量计算的黑匣子。
2. 共线方程与旋转矩阵：外方位元素如何定义一束光的姿态
2.1 共线条件：一次成像的几何本质
任何摄影测量的解算都从共线条件出发：物方点、镜头中心、像点三者位于同一条直线上。用公式写出来就是：
x - x0 = -f * [a1*(X-Xs) + b1*(Y-Ys) + c1*(Z-Zs)] / [a3*(X-Xs) + b3*(Y-Ys) + c3*(Z-Zs)]
y - y0 = -f * [a2*(X-Xs) + b2*(Y-Ys) + c2*(Z-Zs)] / [a3*(X-Xs) + b3*(Y-Ys) + c3*(Z-Zs)]
其中 X、Y、Z 是物方点坐标，Xs、Ys、Zs 是摄影中心坐标，也就是外方位线元素；a1 到 c3 是旋转矩阵的九个元素，由外方位角元素组合而来；f、x0、y0 是内方位元素。这个方程的几何含义非常直接：先把物方点从物方坐标系旋转到像空间坐标系，得到它在相机坐标里的三个分量，然后按针孔投影关系除以 Z 分量并乘以 -f，落到像平面上。
这里有一个经常被忽略的机制细节：旋转矩阵不能随便找九个角度函数拼出来。它必须满足正交性，且行列式为 +1，否则坐标变换会引入尺度和镜像变形。九个元素实际上只有三个自由度，正好对应三个角元素。因此整个平差过程中，矩阵的九个元素不能当作九个独立未知数去修改，必须始终保持由角元素重新生成——这也是为什么每次迭代都要用当前角元素重算旋转矩阵，而不是在矩阵元素上直接加增量，加了就破坏正交性，几何模型立刻失真。
2.2 旋转矩阵：转角系统的定义决定了代码里乘法的顺序
在 VC++ 工程里，计算旋转矩阵几乎是每个平差模块的起点。以摄影测量教科书里最常见的 φ-ω-κ 转角系统为例，它的定义是：先绕 Y 轴转 φ、再绕 X 轴转 ω、最后绕 Z 轴转 κ，旋转矩阵写成：
// 由 phi-omega-kappa 转角系统计算旋转矩阵 R（物方 -> 像空间） int ComputeRotationMatrix(double phi, double omega, double kappa, double R[3][3]) { double cosPhi = cos(phi), sinPhi = sin(phi); double cosOmega = cos(omega), sinOmega = sin(omega); double cosKappa = cos(kappa), sinKappa = sin(kappa); R[0][0] = cosPhi * cosKappa - sinPhi * sinOmega * sinKappa; R[0][1] = -cosPhi * sinKappa - sinPhi * sinOmega * cosKappa; R[0][2] = -sinPhi * cosOmega; R[1][0] = cosOmega * sinKappa; R[1][1] = cosOmega * cosKappa; R[1][2] = -sinOmega; R[2][0] = sinPhi * cosKappa + cosPhi * sinOmega * sinKappa; R[2][1] = -sinPhi * sinKappa + cosPhi * sinOmega * cosKappa; R[2][2] = cosPhi * cosOmega; return 0; }这段代码里最需要留意的是 -sinOmega 出现在 R[1][2] 的位置，这是 Y-X-Z 固定顺序旋转的产物。如果换了别的转角系统，比如先绕 X 轴再绕 Y 轴的 ω-φ-κ 系统，矩阵里的符号分布就会整体变化。我曾经在一个模拟项目X里从某个开源算法里抄来旋转矩阵，却直接套用另一套转角系统输出的角元素，结果平差迭代怎么都不收敛。排查这种问题有很直接的办法：把已知物方坐标的点正算回影像，看算出来的像点是否和观测一致。能对上说明旋转矩阵和角元素定义是配套的；对不上就是转角系统错配，后面所有计算都是白费。
2.3 从物理方程到待估参数：平差要解的未知数到底是什么
光束法把两类未知数同时纳入解算：一类是每张影像的六个外方位元素 Xs、Ys、Zs、φ、ω、κ，它们决定光束在空中的位置和姿态；另一类是每个待定物方点的三个坐标 X、Y、Z。如果影像有 m 张、待定点有 n 个，总未知数就是 6m + 3n。一个像点观测可以提供 x、y 两个方程，因此只要观测总数大于未知数总数，就具备平差条件。
但这里有个入门者容易混淆的点：共线方程本身是非线性的，不能直接列成线性最小二乘。解法是把未知数拆成"近似值加改正数"，围绕当前近似值做泰勒展开，保留一阶项，得到一个关于改正数的线性方程组，然后迭代求解。这个过程称为线性化。每次迭代都重新计算旋转矩阵、重算残差，直到改正数小到阈值以下。也就是说，光束法的计算流程本质上是"线性化、解算、更新、再线性化"的循环，而不是一次性列出方程解完。理解这一点，就理解了为什么平差代码总是带一个迭代主循环，也理解了为什么初值质量能直接决定最终成败。
3. 误差方程与法方程：区域网平差的矩阵骨架
3.1 误差方程：每个像点贡献两行系数
对共线方程线性化之后，每个像点观测可以写成：
vx = a1ΔXs + a2ΔYs + a3ΔZs + a4Δφ + a5Δω + a6Δκ + b1ΔX + b2ΔY + b3*ΔZ - lx
vy 同理。写成矩阵形式就是 V = AΔE + BΔG - L。这里 A 是 2×6 系数矩阵，对应外方位元素改正数；B 是 2×3 系数矩阵，对应物方点坐标改正数；L 是观测值减去当前近似值正算结果的常量项。A 和 B 里每一项本质上都是共线方程对相应未知数的偏导数。
代码实现里最稳妥的做法是解析求导。以 ∂x/∂Xs 为例，令 Xbar、Ybar、Zbar 为物方点坐标经旋转矩阵变换后的三个分量，则 x = -f * Xbar / Zbar，于是：
∂x/∂Xs = f * (R[0][0] * Zbar - R[2][0] * Xbar) / Zbar²
对物方坐标的偏导有个更省事的性质：因为 dX = X - Xs，所以 ∂x/∂X = -∂x/∂Xs，∂x/∂Y = -∂x/∂Ys，∂x/∂Z = -∂x/∂Zs。代码里可以复用已经算好的外方位偏导，把对位置元素的偏导取个负号就能得到对物方坐标的偏导。对三个角元素的偏导要稍微绕一些，需要先求旋转矩阵元素对 φ、ω、κ 的导数，再组合进共线方程。实用的工程技巧是：先完成对 Xs、Ys、Zs 的偏导，再对 X、Y、Z 取相反数，角元素部分单独封装一个函数，这样出错时容易定位，调试时也能逐块验证。
3.2 权阵与最小二乘：不是所有观测都值得同等信任
误差方程建立之后，就要按最小二乘准则加权求解，目标函数是 VᵀPV 取极小。权阵 P 反映每个像点观测的可靠程度。常规做法是以单位权中误差为基准，按像点坐标量测精度给权；控制点或连接点如果有先验精度信息，也可以把它们作为带权约束的附加观测方程加进去。
权阵的选取直接影响平差结果。我有一次把连接点的权设得过高，导致区域网为了迁就少量连接点而把影像姿态整体拉得扭曲，残差看着变小了，接边却出现明显的错位。后来改成按点位实际量测精度估计权，问题才消失。权重不是锦上添花，它是平差模型里一个跟观测流程强相关的自由度。如果影像匹配的同名点精度高，权就大；如果是人工量测或低纹理区域的匹配点，权就要相应下调，否则整体解会被少数噪声大的点带偏。
3.3 法方程的组装与分块：区域网为什么能解上千个未知数
按最小二乘准则，误差方程 V = AΔE + BΔG - L 对应的法方程组装起来是：
| AᵀPA AᵀPB | | ΔE | | AᵀPL | | BᵀPA BᵀPB | | ΔG | = | BᵀPL |
直接对这个大型对称矩阵求逆完全不现实。当影像量和待定点量都到几百时，总未知数几千，稠密矩阵的内存和计算量都会爆炸。光束法能成立的关键在于，这个法方程矩阵有非常规则的块稀疏结构。每个待定点只出现在观测到它的那几张影像对应的行块里，因此 BᵀPB 是一个分块对角矩阵，每个待定点对应一个独立的 3×3 块；AᵀPA 由各影像自身的观测累加而成，同样呈分块对角形态。
利用这个结构，标准做法是先做 Schur 消元，从法方程中把待定点坐标增量 ΔG 消去。先解出用 ΔE 表示的 ΔG，再代回原方程，就得到只含外方位元素的简化法方程：
ΔE = (AᵀPA - AᵀPB (BᵀPB)⁻¹ BᵀPA)⁻¹ (AᵀPL - AᵀPB (BᵀPB)⁻¹ BᵀPL)
由于 BᵀPB 是块对角的，(BᵀPB)⁻¹ 可以逐块求 3×3 逆矩阵，计算量和存储量都很小。简化后的法方程维度是 6m × 6m，只跟影像数量相关。几百张影像时也就是几千维，稠密求解器完全扛得住。这就是光束法区域网平差能在 VC++ 工程里落地的数学基础：不是暴力求解，而是靠结构化的分块消元把问题尺寸压到可处理范围。
4. VC++实现：从数据结构到平差核心的完整链条
4.1 数据结构：矩阵类与观测存储
写平差代码，第一件事是把矩阵运算封装好。VC++ 工程里常用的做法是写一个简单的动态二维数组类，支持对称矩阵存储和高斯消元或 Cholesky 分解。不需要引入庞大的数值库，因为法方程对称正定，Cholesky 分解比通用高斯消元快一倍且数值更稳。观测数据的核心结构可以定义成：
// 像点观测结构体 struct Observation { int imageId; // 影像编号 int pointId; // 待定点编号 double x, y; // 像点坐标（已减去像主点 x0, y0） double weight; // 观测权 };影像的外方位元素用一个数组按 imageId 索引对齐，待定点坐标用另一个数组按 pointId 索引。每次迭代时，遍历所有观测，逐条累加误差方程系数到法方程。这个遍历是整个平差最重的循环，合理的做法是一次循环同时更新法方程的各个分块和右端项，避免为同一个观测重复计算旋转矩阵。内存管理上要特别注意，法方程矩阵在迭代中反复读写，用连续内存的一维数组按行优先存储比 vector 快得多。
4.2 核心代码：误差方程系数与法方程累加
下面给出一段关键代码的骨架，展示单个观测如何生成系数并累加进法方程。偏导数这里用了解析表达式中最核心的平移部分，角元素部分的展开留在完整代码里：
// 计算单个像点观测的误差方程系数 // 输入：物方点坐标近似值、外方位元素近似值、像点观测值 // 输出：a[6] 外方位元素偏导，b[3] 物方坐标偏导，lx/ly 残差常量项 void BuildAndAccumulate( double X, double Y, double Z, double Xs, double Ys, double Zs, double phi, double omega, double kappa, double f, double xObs, double yObs, double weight, double Naa[6][6], double Nab[6][3], double Nbb[3][3], double Ua[6], double Ub[3]) { double R[3][3]; ComputeRotationMatrix(phi, omega, kappa, R); double dX = X - Xs, dY = Y - Ys, dZ = Z - Zs; double Xbar = R[0][0]*dX + R[0][1]*dY + R[0][2]*dZ; double Ybar = R[1][0]*dX + R[1][1]*dY + R[1][2]*dZ; double Zbar = R[2][0]*dX + R[2][1]*dY + R[2][2]*dZ; double xCalc = -f * Xbar / Zbar; double yCalc = -f * Ybar / Zbar; double lx = xObs - xCalc; double ly = yObs - yCalc; // 对平移元素的偏导 double a[6], b[3]; a[0] = f * (R[0][0] * Zbar - R[2][0] * Xbar) / (Zbar * Zbar); // dX/dXs a[1] = f * (R[0][1] * Zbar - R[2][1] * Xbar) / (Zbar * Zbar); // dX/dYs a[2] = f * (R[0][2] * Zbar - R[2][2] * Xbar) / (Zbar * Zbar); // dX/dZs // y 方向的六个偏导、物方坐标反向偏导、角元素偏导在此处补充 // 累加到法方程 for (int i = 0; i < 6; ++i) { Ua[i] += weight * a[i] * lx; for (int j = 0; j <= i; ++j) { Naa[i][j] += weight * a[i] * a[j]; } for (int j = 0; j < 3; ++j) { Nab[i][j] += weight * a[i] * b[j]; } } for (int i = 0; i < 3; ++i) { Ub[i] += weight * b[i] * lx; for (int j = 0; j <= i; ++j) { Nbb[i][j] += weight * b[i] * b[j]; } } }这段代码展示的是单方向的累加逻辑，实际工程要对 x、y 两个方向分别计算系数、分别累加。函数名最开始的 Build 部分是纯计算，不修改外部状态；Accumulate 部分只做对称矩阵的加法。把这两件事分开，是后续调试的法宝——可以单独验证偏导数公式，再验证累加逻辑。我在这里吃过亏：最初把两步混在一个循环里，偏导数公式错了一位，但日志里看到的错误在累加部分，排查了整整一天。拆开之后，用一组已知外方位元素和物方坐标做数值差分验证偏导，十分钟就定位了问题。
4.3 迭代流程：从近似值到收敛
完整的迭代流程可以概括为几步。第一步，用近似外方位元素和近似物方坐标正算所有像点，得到残差向量 L。第二步，按当前近似值计算每个观测的 A、B 系数，累加法方程。第三步，对法方程做 Schur 消元，得到简化法方程，用 Cholesky 分解求外方位元素改正数 ΔE。第四步，回代到分块对角方程，求每个待定点的 ΔG。第五步，用阻尼系数控制步长，更新所有未知数，重新正算残差。重复以上步骤，直到最大改正数小于阈值或者迭代次数达到上限。
迭代控制上有三个参数最容易踩坑。收敛阈值不能只看残差总量，残差总量可能已经被平均掉了，要看最大单个改正数与前一次迭代相比的下降比例。阻尼系数 λ 的初值和调整策略直接决定 LM 迭代是稳定收敛还是振荡发散，初值取 1e-3 并由收敛情况按十倍增减是常用的起点。迭代次数上限设在 50 到 100 次之间，超过这个范围还收敛不了，说明模型本身有问题，靠增加迭代次数没有意义。实际中大多数正常区域网十几次迭代就能收敛，如果三十次还在缓慢下降，应该回头检查初值、连接点分布和权阵设置。
5. 平差避坑指南：初值、病态与发散排查
5.1 发散：迭代次数增加时残差反而越来越大
现象：程序跑起来之后，第一次迭代残差很大，第二、三次更大，甚至直接溢出成 NaN。
原因：最常见的是外方位角元素初值离真值太远。共线方程线性化只在一阶近似有效的邻域内成立，角元素差十几度时，泰勒展开的截断误差已经盖过梯度方向，迭代自然偏离。旋转摄影和倾斜摄影数据里尤其常见，因为影像姿态变化剧烈，初值给零度几乎等于乱猜。
解决：先用单像空间后方交会或用 POS 记录的姿态角做初值，确保所有角元素初值与真值偏差控制在几度以内。同时把每一步迭代的最大改正数打印出来，如果某个角元素改正数超过 0.1 弧度，几乎可以断定是初值问题，停止迭代回到初值检查，而不是继续往下跑。
5.2 病态：法方程条件数过大导致解出来的改正数像噪声
现象：矩阵分解成功，解出来的改正数数值很大且正负交替，看起来完全不符合物理规律。残差可能下降，但外方位元素和真值的差距反而更大。
原因：控制点分布太集中、待定点交会角过小，或者某些影像之间几乎没有连接点，都会让法方程接近奇异。弱约束区域表现为法方程有接近零的特征值，条件数可能到 10 的十二次方以上，数值解自然不稳定。
解决：检查每张影像的连接点数量和空间分布，保证区域网在每块区域都有足够几何约束。必要时对法方程做对角预条件，把六个外方位元素的量纲统一。位置元素是米量级、角元素是弧度量级，直接混在一起装配法方程，数字本身差三个数量级，加预条件或先做单位换算能显著改善数值稳定性。
5.3 旋转矩阵错配：正算对得上，平差却始终不收敛
现象：用已知点正算像点，结果和观测值能对上几毫米；但一进平差迭代就发散，或者外方位元素在合理值附近来回振荡。
原因：绝大多数情况是旋转矩阵与角元素转角系统不匹配。正算能对上是因为测试时用了同一套代码；而平差里的角元素增量方向、偏导数公式用的是另一套转角系统定义，两边对不齐，梯度方向就全错了。
解决：在代码注释里写明转角系统的旋转顺序，并加一个自检函数：输入一组已知外方位元素，正算若干点，再把外方位元素绕单个轴转一小角度，观察像点移动方向是否与偏导数预测方向一致。这个自检能通过，转角系统配准问题就不会在平差阶段暴雷。
5.4 局部极小：残差收敛了，接边还是错位
现象：迭代正常收敛，单位权中误差看起来也不错，但区域网接边区域的同名点错位明显，或者网形出现整体扭曲。
原因：非线性最小二乘只能保证局部收敛，初值落进错误的吸引域时，残差再小也是错误解。常见诱因是连接点数量不足，或者连接点只分布在影像边缘，区域网缺乏内部稳定约束，平差结果就呈现出整片橡皮布被拉抻的变形。
解决：统计每张影像的连接点数量和分布，保证均匀覆盖。对大型区域网，先做一次不带地面控制点的自由网平差，再与已有控制点做相似变换，这套流程能避开不少局部极小问题。如果条件允许，用不同初值跑两遍，比较两组解的物方点坐标差，差异大说明存在多解性。
5.5 精度评估失真：中误差看起来好，实际成果不能用
现象：所有指标都显示平差成功，但最终成果与实地测量对比相差很大，而单位权中误差却小得诱人。
原因：常见的原因有两个。一是把控制点当作无误差的固定值，而不是作为带权观测，控制点自身的量测误差没有进入模型，会扭曲周围网形；二是自由度计算错误，误差方程的观测总数没有减去被消掉的未知数个数，导致单位权中误差系统性偏小。
解决：把控制点作为带权观测加入误差方程，保留它们自身的误差项。平差报告的像点残差 RMS 和单位权中误差要按正确的自由度计算。合格的光束法成果，像点残差 RMS 应该与像点量测精度处于同一数量级；如果十倍优于量测精度，说明自由度算错了或者权阵设置有问题，需要重新检查。
6. 设计取舍：从稠密法方程到现代图优化BA
经典 VC++ 版光束法之所以到几万未知数就吃力，根因在于简化法方程的稠密存储。6m × 6m 的矩阵，几百张影像还能用，到几千张影像做大规模三维重建时，稠密 Cholesky 分解的时间和内存都不可接受。现代 SLAM 与重建领域的 BA 走了不同的工程路线：把相机和物方点建模成二分图，利用稀疏模式做稀疏 Cholesky 分解或多波前分解，再配合 LM 阻尼的信任域策略。这本质上是求解同一组法方程，只是把稀疏性从数学结构贯彻到了数据结构层面。
取舍上值得留意的一点是：对摄影测量这种以控制点、连接点为主，几千路光束的问题，经典分块消元加稠密求解器依然是最稳妥的方案。代码量小、调试容易、数值行为可预期。真正到了几万张影像的工程级空三，才值得引入稀疏迭代求解器。我自己的习惯是：小区域网先跑一遍经典实现，验证所有系数和偏导数的正确性，再考虑性能优化。直接上现代 BA 框架，一旦出问题，你分不清是数据问题还是求解器问题，那才是真正的黑匣子。
另一个常被忽略的取舍是阻尼参数与评价指标的配合。LM 的阻尼作用本质上是把高斯牛顿方向限制在一个可信半径内，它直接决定收敛速度和抗噪声能力，值得单独做一组对比实验来标定。从那以后，我每次做区域网平差，都会强制走一遍转角系统自检、初值检查、连接点分布统计和微型网数值验证，这套流程几乎成了条件反射。面对代码里密密麻麻偏导数时，也只有这些前置检查能给我安全感。希望这些机制层面的拆解，能帮你少踩一些我当年踩过的坑。
本文还有配套的精品资源，点击获取
