做过控制网平差的人早晚会撞上同一个坎手上同时有 GNSS 基线、导线边、水准测段精度明显不在一个档次上全塞进一个平差模型里权该怎么给拍脑袋定一个权比平差报告看着人模人样一到精度评定和粗差检验就露馅——改正数分布明显不对某些观测被系统性地压住。这时候多数人的出路只有一条上 Helmert 方差分量估计让数据自己把每一类观测的单位权方差交出来。Helmert 方差分量估计这件事说穿了很朴素把每类观测各自的单位权方差从已知量降级成未知量和坐标未知数一起解解完再用结果反过来修正权阵来回迭代到收敛。它不挑领域控制测量、变形监测、精密工程测量、甚至实验室多仪器比对的拟合问题都能用但它也不是随手就能上的迭代发散、负方差分量、S 矩阵奇异这几个坑几乎每个第一次上手的人都踩过。下面我按自己做项目的顺序把动机、公式、代码实现和排错流程完整走一遍刚接触平差的新手能看懂思路写过平差程序的老手可以直接抄代码。1. 多源观测混着平差为什么总要栽在权上1.1 一个几乎人人都犯的定权假设标准平差教材里权被定义成 $P_i \sigma_0^2 / \sigma_i^2$其中 $\sigma_0^2$ 是单位权方差$\sigma_i^2$ 是第 $i$ 类观测的方差。公式很干净问题在于 $\sigma_i^2$ 通常根本不知道。多类观测混在一起时最常见的做法是挑一类看起来最靠谱的观测默认它的单位权方差等于 1其他类按经验公式折算权重。水准网就按测段长度定权$P C/S$GNSS 基线要么按边长给个经验公式要么干脆全给等权。这类经验公式背后都藏着一个假设同类观测的误差在统计意义上是同质的而且不同类之间的比例关系可以用一个简单函数描述。假设本身不算离谱但现实里这个假设经常塌方——同一段水准路线柏油路和山间小路每公里的精度可能差一倍同一批 GNSS 基线短基线和长基线、不同解算软件输出的标称精度量级差异能到几倍。1.2 权比错了参数会偏多少这里有个流传很广的误解很多人以为权定错了坐标就会偏。其实不是。只要函数模型没错、观测值本身没粗差最小二乘的参数估值在权阵任意正定的情况下都保持无偏权只影响估计的效率也就是参数的方差大小。真正被权比拖垮的是另外三样东西。一是参数解不再是最优的方差达不到 Cramér–Rao 下界高精度观测本来该有的贡献被稀释二是单位权方差 $\hat{\sigma}_0$ 估计失真进而整份精度评定报告都不可信三是假设检验失效粗差探测用的 $\chi^2$ 统计量和标准化残差前提都是权阵正确权比一错检验的显著性水平就不是你设定的那个了。做过实测的人应该有体会权比不对时粗差定位经常指错点明明是个粗差检验却把旁边那个正常观测标红。这才是最要命的地方。1.3 Helmert 估计到底解决了什么问题Helmert 估计的思路是把上面那个已知的 $\sigma_i^2$变成待估参数。先给一个粗糙的先验权平差一次拿到改正数向量 $V$然后看每一类观测的 $V_i^T P_i V_i$如果这一类观测的权给得偏大它的加权残差平方和就会偏大反之偏小。用这个信号反推出每一类真实的单位权方差再回去修正权阵循环几轮直到各组的单位权方差都落在 1 附近。它和最小二乘的关系可以类比成用平差结果反过来校正平差前提。整个过程不需要额外的外部精度信息靠的是观测数据自身的内部一致性。代价是它对方差分量的可估性有要求样本量太小时估计结果会剧烈抖动这也是后面排错章节要重点讲的内容。理解了这一层动机再去看那些迹运算公式就不会觉得是天书了。2. Helmert 方差分量估计的公式是怎么推出来的2.1 分组误差方程与法方程假设有 $m$ 类观测第 $i$ 类观测的误差方程为$$V_i B_i \hat{x} - l_i, \quad 权阵 P_i, \quad 观测量 n_i$$把所有类堆起来得到整体形式 $V B\hat{x} - l$其中权阵是分块对角的 $P \mathrm{diag}(P_1, P_2, \dots, P_m)$。组成法方程$$N \sum_{i1}^{m} B_i^T P_i B_i \sum_{i1}^{m} N_i, \qquad W \sum_{i1}^{m} B_i^T P_i l_i, \qquad \hat{x} N^{-1} W$$这里 $N_i B_i^T P_i B_i$ 是第 $i$ 类观测对法方程矩阵的贡献。注意一个关键点如果所有观测采用的都是相对权那么整体平差时权的绝对尺度无关紧要参数解不会变变的是单位权方差。Helmert 估计就是在这个自由度上做文章。2.2 二次型期望从 $V_i^T P_i V_i$ 到 $S$ 矩阵方差的估计靠的是加权残差平方和的统计性质。设 $D(l_i) \sigma_i^2 P_i^{-1}$也就是先验权 $P_i$ 对应真实单位权方差 $\sigma_i^2$并且各类之间独立。经过一番推导核心是把 $\hat{x} N^{-1}\sum_j B_j^T P_j l_j$ 代入二次型再逐项求期望可以得到$$\mathrm{E}\left[V_i^T P_i V_i\right] \sigma_i^2\left(n_i - 2,\mathrm{tr}(N^{-1}N_i)\right) \sum_{j1}^{m} \sigma_j^2, \mathrm{tr}(N^{-1}N_i N^{-1}N_j)$$把这个式子按 $i 1, \dots, m$ 排成矩阵形式就得到经典的结构矩阵$$S_{ij} \mathrm{tr}(N^{-1}N_i N^{-1}N_j) \delta_{ij}\left(n_i - 2,\mathrm{tr}(N^{-1}N_i)\right)$$其中 $\delta_{ij}$ 是克罗内克符号。再令 $\theta_i V_i^T P_i V_i$整组方程写成$$S \hat{\sigma}^2 \theta$$解之即得各组的单位权方差估计 $\hat{\sigma}^2 [\hat{\sigma}_1^2, \dots, \hat{\sigma}_m^2]^T$。这套公式有个很漂亮的自我验证如果所有 $\sigma_j^2$ 都等于同一个 $\sigma^2$那么 $\mathrm{E}[V_i^T P_i V_i] \sigma^2(n_i - \mathrm{tr}(N^{-1}N_i))$对 $i$ 求和利用 $\sum_j N_j N$立刻得到 $\mathrm{E}[V^T P V] \sigma^2(n - t)$正是经典的单位权方差公式。所以 Helmert 估计不是另起炉灶它就是经典公式在多类观测下的自然推广。2.3 迭代解算格式与收敛判据$\theta$ 里含的是残差而残差本身依赖于权。所以 $S\hat{\sigma}^2 \theta$ 不能一次解完必须迭代。标准流程是给出一组先验权 $P_i^{(0)}$通常凭经验或标称精度做一次整体平差得到 $\hat{x}$ 和各类残差 $V_i$构造 $S$ 和 $\theta$解出 $\hat{\sigma}_i^2$更新权阵 $P_i^{(k1)} P_i^{(k)} / \hat{\sigma}_i^2$也就是把单位权方差归一化到 1回到第 2 步直到 $\max_i |\hat{\sigma}_i^2 - 1|$ 小于给定阈值常见取 $10^{-6}$ 到 $10^{-8}$。这里有个细节容易搞混更新权时是除以估计值不是乘以。原因很直观——如果第 $i$ 类观测估出来的单位权方差大于 1说明它相对于当前权来说精度比假定的差权应该往下调所以是除以大于 1 的数。收敛后所有类的 $\hat{\sigma}_i^2$ 都趋近 1此时权阵就是自洽的。提示不要用相邻两次 $\hat{\sigma}_i^2$ 差值作为唯一判据。初值特别离谱时前几轮变化幅度天然就大容易误判。同时看是否接近 1和变化是否停滞两个条件更稳。3. 用 Python 从零实现一套 Helmert VCE3.1 数据组织B、l、P 三件套怎么分组代码实现的第一步是数据组织。Helmert 估计天然按组工作所以最自然的结构就是三个列表B_list、l_list、P_list每个列表的下标对应一类观测长度等于类数 $m$。组怎么分取决于你的业务GNSS 网可以按卫星系统分GPS 一组、BDS 一组也可以按基线类型分同步环、异步环水准网按等级分多源融合数据按传感器分。要特别注意分组不是越细越好。分组数和观测数的关系决定了 $S$ 矩阵能不能稳定求逆一般要求每组的观测量至少是未知参数个数的几倍组数控制在 2 到 5 之间最稳。我见过有人把 9 类观测塞进去结果 $S$ 矩阵接近奇异解出来的方差分量全在负数区间乱跳。分组数量是要克制的。3.2 单次平差与 S 矩阵的代码落地核心循环不复杂难的是迹运算容易写错维度。下面这段是可以直接跑的完整实现import numpy as np def helmert_vce(B_list, l_list, P_list, max_iter100, tol1e-8, verboseTrue): m len(B_list) n [li.shape[0] for li in l_list] P [np.array(p, dtypefloat) for p in P_list] t B_list[0].shape[1] history [] sigma2 np.ones(m) for k in range(max_iter): # 1) 组装法方程并求解 N np.zeros((t, t)) W np.zeros(t) Ni_list [] for i in range(m): Ni B_list[i].T P[i] B_list[i] Ni_list.append(Ni) N Ni W B_list[i].T P[i] l_list[i] x np.linalg.solve(N, W) Ninv np.linalg.inv(N) # 2) 构造 theta 与 S theta np.zeros(m) S np.zeros((m, m)) for i in range(m): Vi B_list[i] x - l_list[i] theta[i] float(Vi.T P[i] Vi) for j in range(m): S[i, j] np.trace(Ninv Ni_list[i] Ninv Ni_list[j]) S[i, i] n[i] - 2.0 * np.trace(Ninv Ni_list[i]) # 3) 解方差分量 sigma2_new np.linalg.solve(S, theta) history.append(sigma2_new.copy()) if verbose: print(fiter {k:2d} sigma2 {np.round(sigma2_new, 6)}) # 4) 更新权阵归一化单位权方差 for i in range(m): if sigma2_new[i] 0: P[i] P[i] / sigma2_new[i] if np.max(np.abs(sigma2_new - 1.0)) tol: sigma2 sigma2_new break sigma2 sigma2_new return x, sigma2, P, history有两个地方容易翻车。第一S[i, i]的修正项必须在双层循环外面加如果写进内层j循环里对角线会被重复累加 $m$ 次。第二Ninv只求一次逆就够了别在循环里反复调用np.linalg.inv那是纯粹的浪费。3.3 迭代主循环里的收敛监测上面代码里我特意留了history记录每一轮的方差分量这个小改动在实际调试中非常值钱。因为方差分量的收敛过程经常不是单调的前两轮可能在 0.3 和 3 之间来回跳第三轮才开始收。把history画成折线一眼就能看出来是在收敛还是在震荡。如果是震荡说明分组或初值有问题如果是缓慢逼近那就加大迭代上限慢慢等。另外判断收敛时我用的是 $\max_i |\hat{\sigma}_i^2 - 1|$。有人喜欢用 $\log$ 比值作为判据效果差不多但要注意如果某一轮解出来负值取对数会直接报错所以负值必须单独处理。实际工程里我一般会在解出 $\sigma_2$ 后加一道检查任何一个分量小于 0就不更新权阵直接跳出并报警让用户回去检查数据和分组。3.4 一个可以手算验证的小算例光看代码还是虚找个小例子手推一遍心里就有底了。假设同一个未知量 $x$ 被两组仪器各测了 6 次真值 100.0000第一组精度高单位权标准差约 0.5 mm第二组精度低约 2 mm。观测数据的偏差分别是第一组$0.6, -0.7, 0.5, -0.4, 0.3, -0.3$单位 mm第二组$1.5, -2.0, 2.5, -1.0, 1.8, -2.8$单位 mm初始权故意给错$P_1 I_6$$P_2 0.25 I_6$。因为两组偏差之和都是 0法方程解出 $\hat{x} 100.0000$与权无关。第一轮计算$N_1 6$$N_2 1.5$$N 7.5$$\mathrm{tr}(N^{-1}N_1) 0.8$$\mathrm{tr}(N^{-1}N_2) 0.2$$S_{11} 0.8^2 6 - 2\times0.8 5.04$$S_{22} 0.2^2 6 - 2\times0.2 5.64$$S_{12} 0.16$$\theta_1 1.44$$\theta_2 0.25 \times 24.58 6.145$解方程得 $\hat{\sigma}_1^2 0.2514$$\hat{\sigma}_2^2 1.0824$这两个值对照一下真值第一组真实单位权方差是 $0.25$第二组相对当前权是 $1.0$估得相当准。接下来更新权阵$P_1 \leftarrow I/0.2514 3.978I$$P_2 \leftarrow 0.25I/1.0824 0.2310I$继续迭代轮次$\hat{\sigma}_1^2$$\hat{\sigma}_2^2$$P_1$$P_2$00.25141.08241.0000.250011.13500.95353.9780.231020.99631.00343.5050.242231.00001.00003.5180.2414三轮就收敛了。最终的权比 $P_1 : P_2 \approx 14.6 : 1$跟真值 $16 : 1$ 很接近。这个例子里每组只有 6 个观测样本量偏小还能有这个精度说明只要分组合理、模型正确Helmert 估计的收敛性是可以信赖的。想看更多中间结果的可以把 3.2 节的代码和数据接起来跑history会完整打印出来。4. 工程场景里的落地用法4.1 GNSS 网不同系统、不同基线类型的权比GNSS 基线网的定权是 Helmert 估计最经典的战场。基线解算软件通常会输出每条基线的标称精度 $\sigma a b \cdot d$其中 $d$ 是边长$a$、$b$ 是经验系数。这个公式只反映信号质量层面的精度实际误差里还包含对流层延迟、多路径、天线相位中心偏差等系统性因素不同系统、不同观测时段的数据之间可能存在明显的精度差异。我的做法通常是按系统 时段分两组或三组GPS 一组BDS 一组如果有长时段和短时段混用再按观测时长分。先给等权跑一次看看各组的加权残差平方和比例再用这个比例作为初值去迭代。这样初值不会太离谱一般 4 到 6 轮就收敛。有个容易被忽略的点如果网里有多余观测但图形结构不好比如短边环绕长边$N$ 接近病态$S$ 矩阵的条件数会非常大这时解出的方差分量很不可靠需要先做图形优化或者加基准约束。4.2 水准网按测段长度定权后的重估水准测量的经验定权是 $P C / S$$S$ 是测段长度。这个公式假设每公里水准测量的中误差是常数但在山区和平原混合的线路上这个假设明显站不住脚。用 Helmert 估计可以按地形类型或者测量等级分组把每类的单位权方差反算出来再按 $P^{new} P^{old}/\hat{\sigma}_i^2$ 更新。实际操作里有个技巧水准网的分组不要按测段分要按照测段属性分。比如同样是二等水准山区段和平原段各作为一组而不是每一段都单独估计——每段只有 2 个观测后视、前视根本凑不出能求逆的 $S$ 矩阵。我见过有人按路线分组的结果组数比未知数还多直接算崩。4.3 和其他方差分量估计方法怎么选方差分量估计不止 Helmert 一家工程上常用的还有最小范数二次无偏估计MINQUE、极大似然估计MLE、限制极大似然REML以及贝叶斯方法。它们的适用场景差别挺大我整理成一张表方便对照方法是否要求无偏计算量对初值敏感适用场景Helmert 迭代是中中常规控制网、水准网参数少观测多MINQUE是且方差最小高低理论上限方案小样本更稳MLE / REML否有偏但方差小很高高大规模、需考虑非负约束贝叶斯估计否取决于先验低有可靠先验信息的融合场景Helmert 迭代是 MINQUE 在特定先验下的迭代实现两者在很多场景下结果接近但 MINQUE 不依赖迭代、直接给出解析解代价是必须完整计算所有迹项矩阵维度一大就吃不消。对于常规的几百个观测、几十个参数的网Helmert 迭代的性价比最高这也是它在大地测量领域流传最广的原因。观测规模上千、类别超过 5 类时我会考虑换成 REML 或者约束优化形式避免 $S$ 矩阵病态。5. 踩坑记录与排查清单5.1 迭代震荡与发散最典型的症状是 $\hat{\sigma}_i^2$ 在 0.2 和 5 之间来回跳永远不收敛。成因有三类按出现频率排序初值权比过于离谱、分组本身不可估、以及 $S$ 矩阵条件数过大。初值问题最好解决把初始权改成等权再跑一遍如果等权下就收敛说明是你原来的经验公式太偏了。分组不可估的情况更隐蔽表现为迭代开始时收敛很快但到某一轮突然开始震荡检查一下是不是某几组的 $N_i$ 之间有线性相关。有一招亲测有效给更新后的权阵做一次归一化比如固定 $P_1$ 不变当作基准组只更新其他组。这样可以让 $S$ 矩阵的尺度保持稳定避免权重整体漂移带来的数值震荡。做法是在每轮更新后执行P[i] P[i] * (P[0][0,0] / P_new[0][0,0])这种尺度对齐效率不高但确实治病。5.2 恼人的负方差分量方差分量估计出负值是新手最容易慌的情况。首先要明确一点负的方差估计在数学上是完全可能的因为 $S\hat{\sigma}^2 \theta$ 只是一个线性方程组解的正负由数据和权结构决定并没有强制非负约束。出现负值的常见原因有两个一是某一类观测的数据本身内部一致性太好实际方差远小于先验假设二是观测量太少$S$ 矩阵对角元里的 $n_i - 2\mathrm{tr}(N^{-1}N_i)$ 项接近甚至小于零。处理办法分三步走。第一步检查那一类的观测数如果少于 10 个基本可以判定是样本量问题合并到相邻类别里。第二步检查 $N_i$ 是否和别的组高度相关比如同一台仪器同一时段的观测被拆成两组那它们的方差本来就没法分开估。第三步如果确认数据没问题、只是偶然为负工程上可以直接把它截断到一个小正数比如 $10^{-6}$继续迭代观察后续几轮是否稳定。注意截断只是工程上的权宜之计会破坏无偏性。如果负值反复出现且幅度大说明模型本身有问题别硬扛。5.3 S 矩阵奇异与可估性$S$ 矩阵奇异意味着有些方差分量从数据里根本分不开只能估一个线性组合。最典型的例子是两类观测虽然分开列了但它们对法方程的贡献完全成比例比如同一个网里的两组基线边长和图形结构完全一样只有观测时段不同那这两组的方差分量就不可分别估计。判断方法很简单算一下 $S$ 的条件数超过 $10^{10}$ 就要警惕或者直接看 $N_i$ 之间的相关系数矩阵。还有一个更隐蔽的可估性问题当某些参数的方向上完全没有观测控制时$N$ 本身就奇异$N^{-1}$ 不存在整个公式都失效。这种情况下必须先加基准约束最小约束或符合约束让 $N$ 满秩Helmert 估计才有意义。我踩过一次网里有几个点只有一条基线连着平差时出现秩亏结果方差分量估出来一片乱码折腾了半天才发现是基准问题。5.4 常见问题速查表把上面这些坑整理成一张表出问题时按症状查比从头推导快得多症状最可能原因优先排查项迭代不收敛、来回震荡初值权比太离谱改等权重跑对比history曲线方差分量全为负$S$ 矩阵接近奇异查分组相关性合并高度相关组个别分量为负该类观测太少检查 $n_i$低于 10 考虑合并收敛但值离 1 很远分组不可估或模型缺项检查 $N$ 是否满秩加基准约束迭代轮数特别多初值方向对但尺度偏差大对权阵做尺度归一化平差修正数异常大先验权严重失真用等权先验先做粗差探测再补一个小经验调试阶段把max_iter设为 5先看前几轮的方差分量趋势比一口气跑到 100 轮更容易发现问题。收敛快的方法前 3 轮就能看出苗头要震荡的前 3 轮也已经震起来了。我自己做这类估计的时候有个习惯是先拿一组模拟数据把代码跑通模拟数据里真实方差是已知的能直接验证估计值对不对再到实测数据上跑。实测数据没有真值参照一旦结果异常很难判断是代码错还是数据错。另外Helmert 估计的结果不要当成终点它给出的权比只是一个统计意义上的最优最终的平差报告还要结合图形强度、外业质量记录一起看数据值钱的地方比如高等级控制点附近宁可保守一点给大权。
Helmert方差分量估计实战：多源观测平差定权与Python实现
做过控制网平差的人早晚会撞上同一个坎手上同时有 GNSS 基线、导线边、水准测段精度明显不在一个档次上全塞进一个平差模型里权该怎么给拍脑袋定一个权比平差报告看着人模人样一到精度评定和粗差检验就露馅——改正数分布明显不对某些观测被系统性地压住。这时候多数人的出路只有一条上 Helmert 方差分量估计让数据自己把每一类观测的单位权方差交出来。Helmert 方差分量估计这件事说穿了很朴素把每类观测各自的单位权方差从已知量降级成未知量和坐标未知数一起解解完再用结果反过来修正权阵来回迭代到收敛。它不挑领域控制测量、变形监测、精密工程测量、甚至实验室多仪器比对的拟合问题都能用但它也不是随手就能上的迭代发散、负方差分量、S 矩阵奇异这几个坑几乎每个第一次上手的人都踩过。下面我按自己做项目的顺序把动机、公式、代码实现和排错流程完整走一遍刚接触平差的新手能看懂思路写过平差程序的老手可以直接抄代码。1. 多源观测混着平差为什么总要栽在权上1.1 一个几乎人人都犯的定权假设标准平差教材里权被定义成 $P_i \sigma_0^2 / \sigma_i^2$其中 $\sigma_0^2$ 是单位权方差$\sigma_i^2$ 是第 $i$ 类观测的方差。公式很干净问题在于 $\sigma_i^2$ 通常根本不知道。多类观测混在一起时最常见的做法是挑一类看起来最靠谱的观测默认它的单位权方差等于 1其他类按经验公式折算权重。水准网就按测段长度定权$P C/S$GNSS 基线要么按边长给个经验公式要么干脆全给等权。这类经验公式背后都藏着一个假设同类观测的误差在统计意义上是同质的而且不同类之间的比例关系可以用一个简单函数描述。假设本身不算离谱但现实里这个假设经常塌方——同一段水准路线柏油路和山间小路每公里的精度可能差一倍同一批 GNSS 基线短基线和长基线、不同解算软件输出的标称精度量级差异能到几倍。1.2 权比错了参数会偏多少这里有个流传很广的误解很多人以为权定错了坐标就会偏。其实不是。只要函数模型没错、观测值本身没粗差最小二乘的参数估值在权阵任意正定的情况下都保持无偏权只影响估计的效率也就是参数的方差大小。真正被权比拖垮的是另外三样东西。一是参数解不再是最优的方差达不到 Cramér–Rao 下界高精度观测本来该有的贡献被稀释二是单位权方差 $\hat{\sigma}_0$ 估计失真进而整份精度评定报告都不可信三是假设检验失效粗差探测用的 $\chi^2$ 统计量和标准化残差前提都是权阵正确权比一错检验的显著性水平就不是你设定的那个了。做过实测的人应该有体会权比不对时粗差定位经常指错点明明是个粗差检验却把旁边那个正常观测标红。这才是最要命的地方。1.3 Helmert 估计到底解决了什么问题Helmert 估计的思路是把上面那个已知的 $\sigma_i^2$变成待估参数。先给一个粗糙的先验权平差一次拿到改正数向量 $V$然后看每一类观测的 $V_i^T P_i V_i$如果这一类观测的权给得偏大它的加权残差平方和就会偏大反之偏小。用这个信号反推出每一类真实的单位权方差再回去修正权阵循环几轮直到各组的单位权方差都落在 1 附近。它和最小二乘的关系可以类比成用平差结果反过来校正平差前提。整个过程不需要额外的外部精度信息靠的是观测数据自身的内部一致性。代价是它对方差分量的可估性有要求样本量太小时估计结果会剧烈抖动这也是后面排错章节要重点讲的内容。理解了这一层动机再去看那些迹运算公式就不会觉得是天书了。2. Helmert 方差分量估计的公式是怎么推出来的2.1 分组误差方程与法方程假设有 $m$ 类观测第 $i$ 类观测的误差方程为$$V_i B_i \hat{x} - l_i, \quad 权阵 P_i, \quad 观测量 n_i$$把所有类堆起来得到整体形式 $V B\hat{x} - l$其中权阵是分块对角的 $P \mathrm{diag}(P_1, P_2, \dots, P_m)$。组成法方程$$N \sum_{i1}^{m} B_i^T P_i B_i \sum_{i1}^{m} N_i, \qquad W \sum_{i1}^{m} B_i^T P_i l_i, \qquad \hat{x} N^{-1} W$$这里 $N_i B_i^T P_i B_i$ 是第 $i$ 类观测对法方程矩阵的贡献。注意一个关键点如果所有观测采用的都是相对权那么整体平差时权的绝对尺度无关紧要参数解不会变变的是单位权方差。Helmert 估计就是在这个自由度上做文章。2.2 二次型期望从 $V_i^T P_i V_i$ 到 $S$ 矩阵方差的估计靠的是加权残差平方和的统计性质。设 $D(l_i) \sigma_i^2 P_i^{-1}$也就是先验权 $P_i$ 对应真实单位权方差 $\sigma_i^2$并且各类之间独立。经过一番推导核心是把 $\hat{x} N^{-1}\sum_j B_j^T P_j l_j$ 代入二次型再逐项求期望可以得到$$\mathrm{E}\left[V_i^T P_i V_i\right] \sigma_i^2\left(n_i - 2,\mathrm{tr}(N^{-1}N_i)\right) \sum_{j1}^{m} \sigma_j^2, \mathrm{tr}(N^{-1}N_i N^{-1}N_j)$$把这个式子按 $i 1, \dots, m$ 排成矩阵形式就得到经典的结构矩阵$$S_{ij} \mathrm{tr}(N^{-1}N_i N^{-1}N_j) \delta_{ij}\left(n_i - 2,\mathrm{tr}(N^{-1}N_i)\right)$$其中 $\delta_{ij}$ 是克罗内克符号。再令 $\theta_i V_i^T P_i V_i$整组方程写成$$S \hat{\sigma}^2 \theta$$解之即得各组的单位权方差估计 $\hat{\sigma}^2 [\hat{\sigma}_1^2, \dots, \hat{\sigma}_m^2]^T$。这套公式有个很漂亮的自我验证如果所有 $\sigma_j^2$ 都等于同一个 $\sigma^2$那么 $\mathrm{E}[V_i^T P_i V_i] \sigma^2(n_i - \mathrm{tr}(N^{-1}N_i))$对 $i$ 求和利用 $\sum_j N_j N$立刻得到 $\mathrm{E}[V^T P V] \sigma^2(n - t)$正是经典的单位权方差公式。所以 Helmert 估计不是另起炉灶它就是经典公式在多类观测下的自然推广。2.3 迭代解算格式与收敛判据$\theta$ 里含的是残差而残差本身依赖于权。所以 $S\hat{\sigma}^2 \theta$ 不能一次解完必须迭代。标准流程是给出一组先验权 $P_i^{(0)}$通常凭经验或标称精度做一次整体平差得到 $\hat{x}$ 和各类残差 $V_i$构造 $S$ 和 $\theta$解出 $\hat{\sigma}_i^2$更新权阵 $P_i^{(k1)} P_i^{(k)} / \hat{\sigma}_i^2$也就是把单位权方差归一化到 1回到第 2 步直到 $\max_i |\hat{\sigma}_i^2 - 1|$ 小于给定阈值常见取 $10^{-6}$ 到 $10^{-8}$。这里有个细节容易搞混更新权时是除以估计值不是乘以。原因很直观——如果第 $i$ 类观测估出来的单位权方差大于 1说明它相对于当前权来说精度比假定的差权应该往下调所以是除以大于 1 的数。收敛后所有类的 $\hat{\sigma}_i^2$ 都趋近 1此时权阵就是自洽的。提示不要用相邻两次 $\hat{\sigma}_i^2$ 差值作为唯一判据。初值特别离谱时前几轮变化幅度天然就大容易误判。同时看是否接近 1和变化是否停滞两个条件更稳。3. 用 Python 从零实现一套 Helmert VCE3.1 数据组织B、l、P 三件套怎么分组代码实现的第一步是数据组织。Helmert 估计天然按组工作所以最自然的结构就是三个列表B_list、l_list、P_list每个列表的下标对应一类观测长度等于类数 $m$。组怎么分取决于你的业务GNSS 网可以按卫星系统分GPS 一组、BDS 一组也可以按基线类型分同步环、异步环水准网按等级分多源融合数据按传感器分。要特别注意分组不是越细越好。分组数和观测数的关系决定了 $S$ 矩阵能不能稳定求逆一般要求每组的观测量至少是未知参数个数的几倍组数控制在 2 到 5 之间最稳。我见过有人把 9 类观测塞进去结果 $S$ 矩阵接近奇异解出来的方差分量全在负数区间乱跳。分组数量是要克制的。3.2 单次平差与 S 矩阵的代码落地核心循环不复杂难的是迹运算容易写错维度。下面这段是可以直接跑的完整实现import numpy as np def helmert_vce(B_list, l_list, P_list, max_iter100, tol1e-8, verboseTrue): m len(B_list) n [li.shape[0] for li in l_list] P [np.array(p, dtypefloat) for p in P_list] t B_list[0].shape[1] history [] sigma2 np.ones(m) for k in range(max_iter): # 1) 组装法方程并求解 N np.zeros((t, t)) W np.zeros(t) Ni_list [] for i in range(m): Ni B_list[i].T P[i] B_list[i] Ni_list.append(Ni) N Ni W B_list[i].T P[i] l_list[i] x np.linalg.solve(N, W) Ninv np.linalg.inv(N) # 2) 构造 theta 与 S theta np.zeros(m) S np.zeros((m, m)) for i in range(m): Vi B_list[i] x - l_list[i] theta[i] float(Vi.T P[i] Vi) for j in range(m): S[i, j] np.trace(Ninv Ni_list[i] Ninv Ni_list[j]) S[i, i] n[i] - 2.0 * np.trace(Ninv Ni_list[i]) # 3) 解方差分量 sigma2_new np.linalg.solve(S, theta) history.append(sigma2_new.copy()) if verbose: print(fiter {k:2d} sigma2 {np.round(sigma2_new, 6)}) # 4) 更新权阵归一化单位权方差 for i in range(m): if sigma2_new[i] 0: P[i] P[i] / sigma2_new[i] if np.max(np.abs(sigma2_new - 1.0)) tol: sigma2 sigma2_new break sigma2 sigma2_new return x, sigma2, P, history有两个地方容易翻车。第一S[i, i]的修正项必须在双层循环外面加如果写进内层j循环里对角线会被重复累加 $m$ 次。第二Ninv只求一次逆就够了别在循环里反复调用np.linalg.inv那是纯粹的浪费。3.3 迭代主循环里的收敛监测上面代码里我特意留了history记录每一轮的方差分量这个小改动在实际调试中非常值钱。因为方差分量的收敛过程经常不是单调的前两轮可能在 0.3 和 3 之间来回跳第三轮才开始收。把history画成折线一眼就能看出来是在收敛还是在震荡。如果是震荡说明分组或初值有问题如果是缓慢逼近那就加大迭代上限慢慢等。另外判断收敛时我用的是 $\max_i |\hat{\sigma}_i^2 - 1|$。有人喜欢用 $\log$ 比值作为判据效果差不多但要注意如果某一轮解出来负值取对数会直接报错所以负值必须单独处理。实际工程里我一般会在解出 $\sigma_2$ 后加一道检查任何一个分量小于 0就不更新权阵直接跳出并报警让用户回去检查数据和分组。3.4 一个可以手算验证的小算例光看代码还是虚找个小例子手推一遍心里就有底了。假设同一个未知量 $x$ 被两组仪器各测了 6 次真值 100.0000第一组精度高单位权标准差约 0.5 mm第二组精度低约 2 mm。观测数据的偏差分别是第一组$0.6, -0.7, 0.5, -0.4, 0.3, -0.3$单位 mm第二组$1.5, -2.0, 2.5, -1.0, 1.8, -2.8$单位 mm初始权故意给错$P_1 I_6$$P_2 0.25 I_6$。因为两组偏差之和都是 0法方程解出 $\hat{x} 100.0000$与权无关。第一轮计算$N_1 6$$N_2 1.5$$N 7.5$$\mathrm{tr}(N^{-1}N_1) 0.8$$\mathrm{tr}(N^{-1}N_2) 0.2$$S_{11} 0.8^2 6 - 2\times0.8 5.04$$S_{22} 0.2^2 6 - 2\times0.2 5.64$$S_{12} 0.16$$\theta_1 1.44$$\theta_2 0.25 \times 24.58 6.145$解方程得 $\hat{\sigma}_1^2 0.2514$$\hat{\sigma}_2^2 1.0824$这两个值对照一下真值第一组真实单位权方差是 $0.25$第二组相对当前权是 $1.0$估得相当准。接下来更新权阵$P_1 \leftarrow I/0.2514 3.978I$$P_2 \leftarrow 0.25I/1.0824 0.2310I$继续迭代轮次$\hat{\sigma}_1^2$$\hat{\sigma}_2^2$$P_1$$P_2$00.25141.08241.0000.250011.13500.95353.9780.231020.99631.00343.5050.242231.00001.00003.5180.2414三轮就收敛了。最终的权比 $P_1 : P_2 \approx 14.6 : 1$跟真值 $16 : 1$ 很接近。这个例子里每组只有 6 个观测样本量偏小还能有这个精度说明只要分组合理、模型正确Helmert 估计的收敛性是可以信赖的。想看更多中间结果的可以把 3.2 节的代码和数据接起来跑history会完整打印出来。4. 工程场景里的落地用法4.1 GNSS 网不同系统、不同基线类型的权比GNSS 基线网的定权是 Helmert 估计最经典的战场。基线解算软件通常会输出每条基线的标称精度 $\sigma a b \cdot d$其中 $d$ 是边长$a$、$b$ 是经验系数。这个公式只反映信号质量层面的精度实际误差里还包含对流层延迟、多路径、天线相位中心偏差等系统性因素不同系统、不同观测时段的数据之间可能存在明显的精度差异。我的做法通常是按系统 时段分两组或三组GPS 一组BDS 一组如果有长时段和短时段混用再按观测时长分。先给等权跑一次看看各组的加权残差平方和比例再用这个比例作为初值去迭代。这样初值不会太离谱一般 4 到 6 轮就收敛。有个容易被忽略的点如果网里有多余观测但图形结构不好比如短边环绕长边$N$ 接近病态$S$ 矩阵的条件数会非常大这时解出的方差分量很不可靠需要先做图形优化或者加基准约束。4.2 水准网按测段长度定权后的重估水准测量的经验定权是 $P C / S$$S$ 是测段长度。这个公式假设每公里水准测量的中误差是常数但在山区和平原混合的线路上这个假设明显站不住脚。用 Helmert 估计可以按地形类型或者测量等级分组把每类的单位权方差反算出来再按 $P^{new} P^{old}/\hat{\sigma}_i^2$ 更新。实际操作里有个技巧水准网的分组不要按测段分要按照测段属性分。比如同样是二等水准山区段和平原段各作为一组而不是每一段都单独估计——每段只有 2 个观测后视、前视根本凑不出能求逆的 $S$ 矩阵。我见过有人按路线分组的结果组数比未知数还多直接算崩。4.3 和其他方差分量估计方法怎么选方差分量估计不止 Helmert 一家工程上常用的还有最小范数二次无偏估计MINQUE、极大似然估计MLE、限制极大似然REML以及贝叶斯方法。它们的适用场景差别挺大我整理成一张表方便对照方法是否要求无偏计算量对初值敏感适用场景Helmert 迭代是中中常规控制网、水准网参数少观测多MINQUE是且方差最小高低理论上限方案小样本更稳MLE / REML否有偏但方差小很高高大规模、需考虑非负约束贝叶斯估计否取决于先验低有可靠先验信息的融合场景Helmert 迭代是 MINQUE 在特定先验下的迭代实现两者在很多场景下结果接近但 MINQUE 不依赖迭代、直接给出解析解代价是必须完整计算所有迹项矩阵维度一大就吃不消。对于常规的几百个观测、几十个参数的网Helmert 迭代的性价比最高这也是它在大地测量领域流传最广的原因。观测规模上千、类别超过 5 类时我会考虑换成 REML 或者约束优化形式避免 $S$ 矩阵病态。5. 踩坑记录与排查清单5.1 迭代震荡与发散最典型的症状是 $\hat{\sigma}_i^2$ 在 0.2 和 5 之间来回跳永远不收敛。成因有三类按出现频率排序初值权比过于离谱、分组本身不可估、以及 $S$ 矩阵条件数过大。初值问题最好解决把初始权改成等权再跑一遍如果等权下就收敛说明是你原来的经验公式太偏了。分组不可估的情况更隐蔽表现为迭代开始时收敛很快但到某一轮突然开始震荡检查一下是不是某几组的 $N_i$ 之间有线性相关。有一招亲测有效给更新后的权阵做一次归一化比如固定 $P_1$ 不变当作基准组只更新其他组。这样可以让 $S$ 矩阵的尺度保持稳定避免权重整体漂移带来的数值震荡。做法是在每轮更新后执行P[i] P[i] * (P[0][0,0] / P_new[0][0,0])这种尺度对齐效率不高但确实治病。5.2 恼人的负方差分量方差分量估计出负值是新手最容易慌的情况。首先要明确一点负的方差估计在数学上是完全可能的因为 $S\hat{\sigma}^2 \theta$ 只是一个线性方程组解的正负由数据和权结构决定并没有强制非负约束。出现负值的常见原因有两个一是某一类观测的数据本身内部一致性太好实际方差远小于先验假设二是观测量太少$S$ 矩阵对角元里的 $n_i - 2\mathrm{tr}(N^{-1}N_i)$ 项接近甚至小于零。处理办法分三步走。第一步检查那一类的观测数如果少于 10 个基本可以判定是样本量问题合并到相邻类别里。第二步检查 $N_i$ 是否和别的组高度相关比如同一台仪器同一时段的观测被拆成两组那它们的方差本来就没法分开估。第三步如果确认数据没问题、只是偶然为负工程上可以直接把它截断到一个小正数比如 $10^{-6}$继续迭代观察后续几轮是否稳定。注意截断只是工程上的权宜之计会破坏无偏性。如果负值反复出现且幅度大说明模型本身有问题别硬扛。5.3 S 矩阵奇异与可估性$S$ 矩阵奇异意味着有些方差分量从数据里根本分不开只能估一个线性组合。最典型的例子是两类观测虽然分开列了但它们对法方程的贡献完全成比例比如同一个网里的两组基线边长和图形结构完全一样只有观测时段不同那这两组的方差分量就不可分别估计。判断方法很简单算一下 $S$ 的条件数超过 $10^{10}$ 就要警惕或者直接看 $N_i$ 之间的相关系数矩阵。还有一个更隐蔽的可估性问题当某些参数的方向上完全没有观测控制时$N$ 本身就奇异$N^{-1}$ 不存在整个公式都失效。这种情况下必须先加基准约束最小约束或符合约束让 $N$ 满秩Helmert 估计才有意义。我踩过一次网里有几个点只有一条基线连着平差时出现秩亏结果方差分量估出来一片乱码折腾了半天才发现是基准问题。5.4 常见问题速查表把上面这些坑整理成一张表出问题时按症状查比从头推导快得多症状最可能原因优先排查项迭代不收敛、来回震荡初值权比太离谱改等权重跑对比history曲线方差分量全为负$S$ 矩阵接近奇异查分组相关性合并高度相关组个别分量为负该类观测太少检查 $n_i$低于 10 考虑合并收敛但值离 1 很远分组不可估或模型缺项检查 $N$ 是否满秩加基准约束迭代轮数特别多初值方向对但尺度偏差大对权阵做尺度归一化平差修正数异常大先验权严重失真用等权先验先做粗差探测再补一个小经验调试阶段把max_iter设为 5先看前几轮的方差分量趋势比一口气跑到 100 轮更容易发现问题。收敛快的方法前 3 轮就能看出苗头要震荡的前 3 轮也已经震起来了。我自己做这类估计的时候有个习惯是先拿一组模拟数据把代码跑通模拟数据里真实方差是已知的能直接验证估计值对不对再到实测数据上跑。实测数据没有真值参照一旦结果异常很难判断是代码错还是数据错。另外Helmert 估计的结果不要当成终点它给出的权比只是一个统计意义上的最优最终的平差报告还要结合图形强度、外业质量记录一起看数据值钱的地方比如高等级控制点附近宁可保守一点给大权。
