企业官网建设流程全解析

简介：这是一份面向摄影测量与计算机视觉学习者的光束法区域网平差VC++工程源码，重点演示如何通过最小二乘与束调整解算像片外方位元素，适用于多视图三维重建、无人机测绘等场景的算法入门与实践。压缩包共32个文件，以h头文件和cpp源码为主体，附带原始观测数据、结果数据及Visual C++工程配置文件，总大小约55KB，结构精简、便于直接编译学习。目前已有1320人学习下载。资源中不仅包含BundleAdjustment核心模块与封装好的Matrix矩阵运算类，还提供了光束法数据、resultDegree与resultRad等示例文本，能够帮助读者对照输入输出理解平差迭代过程，掌握旋转角与平移量的解算思路，适合具备一定C++基础、希望深入摄影测量平差原理的开发者参考。

1. 光束法区域网平差：为什么“平差”不是“凑数”

很多人第一次接触摄影测量里的光束法区域网平差，都会觉得它不过是个“用最小二乘把参数凑出来”的过程：反正方程多了，解个超定方程组而已。但真正手写过VC++实现后你会发现，这个说法掩盖了最关键的地方——光束法平差不是“凑数”，而是用共线方程把影像的外方位元素和地面点坐标捆进一个互相牵制的模型里，然后利用稀疏结构高效求解。它难的不是最小二乘本身，而是“几百张影像、几万个点、几十万次观测”同时参与时，如何让矩阵可控、让运算能收敛、让程序不崩。搞懂它，才能真正看懂航测软件里那些“平差失败”和“像点残差超限”的报错到底在说什么。

2. 光束法平差的数学模型：共线方程与外方位元素

2.1 共线方程的物理含义：像点、控制点和投影中心三点一线

光束法平差的“光束”二字，指的是摄影瞬间物方点、镜头投影中心和像点在同一条直线上。这个几何关系是摄影测量的第一性原理，数学表达就是共线方程。设像点坐标为 x、y（已经加过主点偏移和畸变改正），对应地面点坐标为 X、Y、Z，投影中心坐标为 Xs、Ys、Zs，那么理想情况下满足：

x - x0 = -f * [a1*(X-Xs) + b1*(Y-Ys) + c1*(Z-Zs)] / [a3*(X-Xs) + b3*(Y-Ys) + c3*(Z-Zs)]

这里 a1...c3 是由外方位角元素（通常是 ω、φ、κ）组成的旋转矩阵 R 的九个方向余弦。f 是主距，x0、y0 是像主点坐标。写成方程之后你会发现，像点坐标是地面点坐标和外方位元素的非线性函数，所以平差必须迭代线性化，这是整个算法的起点。

在实现里，旋转矩阵很少直接存九个浮点数再参与运算，而是用一个 3×3 数组，由三个角的三角函数组合而成。常见写法是先按 ω、φ、κ 的转序计算，再累乘。这里有个最容易翻车的地方：摄影测量里旋转矩阵的排列约定和计算机图形学不一样，方向余弦 a1、a2、a3 的符号很容易错。某跨平台项目里在自检时发现所有残差都带同一个偏移，最后才定位到 κ 角的正负号搞反了。

2.2 平差变量：外方位元素、内方位元素与地面点坐标如何分层

光束法平差里，未知数按来源分成三类：每张影像的外方位元素（Xs、Ys、Zs、ω、φ、κ），每个地面点坐标（X、Y、Z），以及可能参与自检校的内方位元素（x0、y0、f 和畸变系数）。所以在区域网里，如果影像数 m、地面点数 n，未知数总量大约是 6m + 3n 甚至更多。这个“6m 加 3n”的结构决定了后面矩阵的组织方式：外方位元素是影像编号索引的，地面点是点编号索引的，而像点观测则是“哪张影像、哪个点”的二元组合。

平差时，观测值是像点坐标量测值，而不是“地面点坐标”。很多人写代码时下意识把像点当自变量，把地面点当因变量，这是理解上的一步之遥，却会导致法方程结构完全错乱。正确思路是：每个像点观测都同时牵连一个外方位元素组和一个地面点坐标组，所以误差方程展开后有六项关于影像未知数的偏导，三项关于地面点的偏导，剩下的才是常数项。

2.3 误差方程与权阵：怎么把“观测值”变成线性方程组

共线方程非线性，所以先取近似值做泰勒展开，得到误差方程：

V = A * δt + B * δX - L

其中 V 是残差向量，δt 是外方位元素改正数，δX 是地面点改正数，A 和 B 分别是误差方程对这些参数的偏导数矩阵，L 是像点观测值减去用近似值计算出的像点坐标（也就是当前模型下的“观测残差”）。这里的每一个观测还需要乘以权重，权重矩阵 P 通常由像点坐标量测精度决定，比如 σx = σy = 2 微米（像元大小的一半），那么权就是 1/σ^2。区域网平差中不同影像可能摄影条件不同，单位权中误差也会不一样，所以在代码里权阵常写成对角线型，按观测来源乘一个比例因子。

这个模型看上去简单，但真正程序化时，要特别留意 A、B 矩阵里偏导数的计算顺序。VC++ 实现里很容易在循环里反复调用三角函数，性能极差；常见的优化是每张影像固化旋转矩阵，用常量数组存放，偏导数公式也拆成临时变量复用，避免对同一个 sin、cos 重复求值。

3. 稀疏矩阵与法方程：区域网平差的性能命门

3.1 法方程为什么天生稀疏：像点观测的局部性

把误差方程代入最小二乘后，得到法方程：

(A^T P A) δ = A^T P L

这里的系数矩阵稀疏度惊人。对于任意一张影像上的一个像点，它只连接了“这一张影像的外方位元素”和“这一个地面点坐标”，不会去连接其他影像或者其他地面点。因此，如果我们按“先影像未知数，后地面点未知数”的块结构排布，法方程矩阵会呈现明显的块状稀疏模式：主对角线附近是每张影像和每个点的自相关块，非对角线只有在“影像-点对”存在观测时才出现元素。

这意味着不能用一个 m+n 阶的稠密二维数组去存储，否则几百张影像加几万点的规模直接超出内存。常见做法是只存非零块：外方位元素块 6×6，地面点块 3×3，连接块 6×3 和 3×6。存储结构用三个动态数组分别记录行号、列号和值，或者用 CRS（Compressed Row Storage）稀疏格式。稀疏格式的实现是光束法平差代码的核心资产，因为后面所有性能优化基本都围绕它展开。

3.2 先消地面点还是先消外方位元素：约化相机矩阵

直接解整个法方程也可以，但在摄影测量里更常见的是先消去地面点坐标。理由很简单：区域网平差中地面点数量往往达几万甚至几十万，而影像数量不过几百，未知数收缩后能极大降低后续分解代价。把法方程按“影像块、交叉块、点块”分块后，利用 Schur 补消去点块，得到只包含影像外方位元素未知数的“约化法方程”，这个矩阵维度只有 6m×6m，虽然不像原矩阵那么稀疏，但规模小得多。解出影像未知数后，再回代求出地面点坐标。

这个“先消点、后解相机”的顺序不是拍脑袋定的，它直接对应实际程序流程：先对整个法方程做一次分块线性处理，把点相关的下三角块消隐，再用 Cholesky 分解解约化矩阵。VC++ 代码里通常用两个函数分别完成这个分块消元和回代。如果影像数也到了一千以上，约化矩阵仍然太大，那就得走迭代法（比如共轭梯度）而不是直接分解。

3.3 VC++实现中的稀疏矩阵：指针、索引和局部性

下面是一段典型的稀疏块法方程累加逻辑，展示了如何在循环里把单个像点观测的贡献累加到法方程对应的块位置。

// 观测索引：imgIdx 影像编号，ptIdx 地面点编号 // A 为 2x6 误差方程关于外方位元素偏导，B 为 2x3 关于地面点偏导 // P 为 2x2 权阵（通常为对角） for (int obs = 0; obs < numObs; ++obs) { int i = imgIdx[obs]; // 影像索引 int j = ptIdx[obs]; // 点索引 // 累加影像自块：A^T * P * A for (int r = 0; r < 6; ++r) for (int c = r; c < 6; ++c) { double sum = 0.0; for (int k = 0; k < 2; ++k) sum += A[k][r] * P[k][k] * A[k][c]; N_img[i][r][c] += sum; } // 累加交叉块：A^T * P * B for (int r = 0; r < 6; ++r) for (int c = 0; c < 3; ++c) { double sum = 0.0; for (int k = 0; k < 2; ++k) sum += A[k][r] * P[k][k] * B[k][c]; N_cross[i][j][r][c] += sum; } // 点块、点-影像对称块类似，此处省略 }

这里的 N_img 和 N_cross 都是按索引预先分配好的三维数组还是稀疏容器，取决于测区规模。如果点数是十万级，N_cross 就不能用三维数组整体申请，因为绝大多数“影像-点”组合没有观测，必须用哈希表或排序后的索引数组去存储。初学者常犯的错是为所有可能组合预先开一个矩阵，结果内存翻车。

另一个值得注意的参数是偏导矩阵的排列顺序。上面代码里 A 是 2×6，对应像点坐标的两个方程分别对六个外方位元素求导。有些实现把行数写成 2、列数是 6，有些反过来；法方程累加逻辑里行列顺序也必须同步调整，否则对称性一错，后面的 Cholesky 分解结果全是错误的。

4. 平差迭代：从初值到收敛

4.1 线性化与法化：迭代循环里到底发生了什么

光束法平差不是解一次线性方程就结束。因为共线方程非线性，解出来的改正数只在线性化近似点附近有效，需要把求得的改正数加到近似值上，重新计算偏导数、重新法化、再解，直到改正数小到可忽略。每一次迭代的主要步骤如下：

1. 用当前外方位元素和地面点坐标计算每个像点的预测坐标，得到观测残差 L。
2. 在当前近似值处计算偏导矩阵 A、B，形成误差方程。
3. 累加法方程系数矩阵和常数项。
4. 消去地面点未知数，解出外方位元素改正数。
5. 回代求地面点改正数，更新全部未知数。
6. 判断收敛，不满足则回到步骤 1。

这个过程在 VC++ 中通常用一个外层 while 循环包住，循环条件写为迭代次数上限和改正数阈值。注意：每次迭代都必须用“更新后的参数值”重新计算旋转矩阵 R 和像点预测值，而不是沿用上一轮的偏导数。很多程序第一次迭代正常，第二次直接发散，原因就是把 R 缓存在函数外面，忘了更新。

在代码实现里，单位权中误差也是一个每次迭代都会变的值。单位权中误差的定义是：

sigma0 = sqrt( (V^T P V) / (n - u) )

其中 n 是观测方程个数，u 是未知数个数。这个自由度差值就是“多余观测数”，它决定了平差有没有意义。如果所有方程加起来恰好等于未知数个数，分母为零，单位权中误差无定义，平差退化。

4.2 迭代终止条件：改正数、单位权中误差和收敛阈值

判断收敛经常用两个指标同时卡：最大改正数绝对值小于某个相对阈值，并且这回的单位权中误差变化量已经很小。代码里长这样：

bool converged = false; double prevSigma0 = 0.0; for (int iter = 0; iter < maxIter; ++iter) { // ... 法方程求解，得到 dx 向量 if (norm_dx < 1e-6 * norm_param) { converged = true; } double sigma0 = computeSigma0(V, numObs, numUnknowns); if (fabs(sigma0 - prevSigma0) < 1e-4 && iter > 0) { converged = true; } prevSigma0 = sigma0; if (converged) break; }

这里的 norm_dx 通常取所有改正数平方和的平方根（2-范数），norm_param 取对应参数向量的 2-范数，用相对值判断能避免不同测区尺度差异导致阈值不合适。单位权中误差变化阈值建议取 1e-4 到 1e-6，太小会多跑几次费时间，太大会提前终止导致精度不足。我在实际项目里体会很深的一点是：判断收敛前至少跑两次迭代，因为第一次的改正数大得离谱，直接套小阈值会认为不收敛而继续跑到迭代上限；正确做法是先跑满几次再启用收敛判断。

4.3 粗差检测：标准化残差里藏着错误匹配

光束法平差的精准度高度依赖观测质量。一个像点坐标如果有几微米的粗差，往往会被平差“分摊”到周围所有参数上，导致局部区域扭曲。所以平差过程里必须有粗差检测步骤。最简单的就是标准残差法：计算每个像点的残差 v，除以它的中误差估值为 σv，得到一个标准化残差；如果某个像点的标准化残差绝对值大于 3，就怀疑有粗差，剔除后重新平差。

在迭代里这会形成一个外层循环：先用所有观测平差，计算每个像点的标准化残差，剔除最大超限点，再重新从步骤 1 迭代。注意每次剔除后自由度和法方程结构都会变化，必须重新法化，不能只删掉法方程里的某一行然后继续解，因为残差本身已经因为“污染”而不可靠。这里常见的血泪经验是：粗差检测的阈值设到 2.5 还是 3.0 取决于观测冗余度，冗余高时阈值可以放严一点。

5. 光束法平差常见问题与排查

5.1 初值不准导致迭代发散

现象：第一次迭代后改正数巨大，第二次直接跑出极大值或 NaN。

原因：共线方程线性化只在小范围内成立。如果旋转角初值差了几度或几十度，偏导数会严重失真，迭代不会逐渐逼近，只会越偏越远。

解决：先单独用空间后方交会求每张影像外方位元素的初值，再用前方交会求地面点初值。区域网平差里的“初值计算”不是可选项，而是必选项。还有一个常见补救是阻尼最小二乘（Levenberg-Marquardt），在法方程对角线加一个逐渐减小的阻尼项，但本质还是要有合理的初值。

5.2 控制点不足导致法方程秩亏

现象：解出来的外方位元素在测区边缘漂移，或者迭代无法收敛，程序报“矩阵奇异”或“分解失败”。

原因：区域网平差必须有足够的地面控制点或像控点来确定空间基准。如果一条航带里所有影像都没有控制点，法方程在有观测的点上仍然不满秩，特别是矩阵的零空间对应测区整体的平移、旋转和缩放自由度。

解决：检查控制点数量是否至少满足每一张参与平差的影像有足够约束。更严谨的办法是查看法方程矩阵 Cholesky 分解时对角线元素是否为 0 或接近 0，从而定位哪些未知数不可估。加控制点，或者引入连接点间的高程约束，是常用手段。

5.3 内存溢出：矩阵分配不规范

现象：程序在点数上万后突然崩溃，报内存不足。

原因：用稠密二维数组存储整个法方程，或为所有“影像-点”组合分配交叉块。100 张影像、10 万个点的测区，如果交叉块全开，轻轻松松几个 GB 起步，而实际观测数远小于这个规模。

解决：改用稀疏存储，只保存有观测的结构块。VC++ 里可以先用哈希表收集非零块的索引，再转换成紧密排列的数组，减少内存碎片。打开编译器优化选项，并避免在循环内使用 new/delete 频繁分配。

5.4 坐标系不一致：最容易忽视的“隐性错误”

现象：平差结果残差看起来不错，但最后生成的三维点坐标和已知控制点相差很大。

原因：影像的外方位元素初值可能基于某个局部坐标系，地面控制点用的是另一个椭球或投影系；或者旋转角系统不同（比如东偏北和北偏东定义互换）。

解决：写代码时把坐标系作为参数显式传入，并且在每轮迭代前做一次坐标预检。可以设计一个自检函数：用已知控制点重投影到影像上，如果重投影误差很大，说明坐标系转换还没对上。这一步看起来是“玄学”，其实是最常见的实际翻车点。

5.5 排查手法：盯住残差和改正数

现象：程序能跑完平差，结果说不上不对，但精度验收不过。

原因：没有中间可观测指标，导致问题被藏在最后的结果里。

解决：每次迭代输出三类信息：最大残差向量、单位权中误差、外方位元素改正数范数。如果残差逐渐减小但减得极慢，说明有系统误差或存在粗差；如果改正数符号在正负间振荡，大概率是矩阵顺序或导数符号有错。这个排查步骤写完就放进日常流程，能省大量时间。

6. 设计取舍：手写 VC++ 平差的边界在哪里

光束法区网平差的实现方式有很多，有人用现成的数值库，有人用 Python 脚本快速验证，也有人坚持用 VC++ 从底层手写稀疏矩阵。我自己的体会是：这个领域里“性能”和“可控性”是一对永远需要权衡的矛盾。通用数学库能保证数值算法的正确性，却未必了解摄影测量里“先消点、后消影像”这类特殊顺序；而手写代码虽然调试成本高，却能在面对大测区时精确掌控内存布局和迭代策略。

对个人项目或小测区，VC++ 配合 Eigen 或者类似模板库是现实选择；但对于需要处理连续测区拼接、上千张影像的生产系统，稀疏数据结构、约化相机矩阵和 Schur 补的实现仍然值得自己写一遍。因为只有手写过一次，你才会理解为什么程序在测区边缘总是不稳，为什么控制点一少矩阵就奇异，为什么残差分布会呈现条带状。从那以后，我每次写完平差模块都会强制自己走一遍“检查秩亏、打印残差、验证收敛”三步，而不是直接看结果文件。希望帮到你。

本文还有配套的精品资源，点击获取