做GNSS数据处理的同学十有八九都遇到过这个头疼问题多系统融合解算时GPS、BDS、Galileo的观测值到底该怎么定权伪距和载波相位精度差着几个量级硬塞进同一个法方程里权比给不对解出来的坐标可能偏得离谱。早些年我习惯按高度角模型给个先验权但项目做多了就发现先验模型只能解决一部分问题——不同系统、不同接收机、不同观测类型的实际噪声水平差别很大单靠经验公式很难拿准。后来系统学习了赫尔默特方差分量估计Helmert Variance Component Estimation才算是把“定权”这块真正补上了。这篇文章就把赫尔默特方差分量估计的原理和Python实现完整拆开讲一遍适合刚接触GNSS数据处理的学生也适合已经在做多源融合定位但觉得定权不够严谨的工程师。1. 为什么GNSS数据处理离不开赫尔默特方差分量估计1.1 平差模型的最后一环随机模型定权任何一个最小二乘平差问题都由两部分组成函数模型和随机模型。函数模型描述观测值和未知参数之间的数学关系比如伪距观测方程、载波相位观测方程、双差观测方程这部分大家都很熟。真正的分歧往往出在随机模型上——观测值的方差协方差矩阵怎么给也就是权阵怎么定。权阵的本质是方差协方差矩阵的逆。如果给某一类观测值的权过大相当于告诉平差系统“我非常相信这些观测”那么解算结果就会倾向于迁就这类观测反之权过小则会被“无视”。在GNSS数据处理里最常见的定权方式是高度角定权卫星高度角越高信号传播路径越短、大气延迟误差越小就给越大的权。这个方法简单有效但它只刻画了“卫星几何条件”这个维度无法反映系统间偏差、接收机通道差异、信号频段差异带来的噪声水平变化。赫尔默特方差分量估计解决的就是这个问题不靠经验拍脑袋而是利用最小二乘平差之后的残差二次型反过来估计每一类观测值的方差分量再用估计出的方差分量修正权阵重新平差迭代若干轮后让权阵与实际观测噪声水平相匹配。这个过程可以类比成团队协作先按某个初步判断分配任务权重干完一阶段看看每组成员的成果质量再重新调整下一阶段的任务分配几轮下来权重就合理了。1.2 赫尔默特估计最典型的几个应用场景赫尔默特方差分量估计在GNSS领域最常见的应用场景是以下几类理解了场景就理解了它的工程价值。第一个就是多系统GNSS融合定位。GPS、BDS、Galileo、GLONASS四类系统观测值联合平差时不同系统的码伪距和载波相位噪声水平并不一致直接用等权或者统一的高度角模型处理往往会造成“强系统带偏弱系统”的问题。用赫尔默特估计可以自动估计各系统观测值之间的权比。第二个是GNSS/INS组合导航。组合导航的观测方程里既有GNSS位置约束又有惯性导航系统的预测值两类观测量的量纲和噪声特性完全不同。惯导在短时间内精度很高长时间会漂移这种时变特性使得固定权比很难适配用方差分量估计做自适应定权是工程上比较成熟的做法。第三个是传统大地测量中的边角网平差。测角和测边的精度单位完全不同角度观测的方差和距离观测的方差不在一个量纲上经典做法就是赫尔默特方差分量估计教科书里的例子几乎都来源于此。第四个是变形监测中的多传感器数据融合。全站仪、GNSS接收机、静力水准仪各自的噪声特性差异很大联合平差时也需要方差分量估计来平衡各类传感器的贡献。总的来说只要你的平差模型里有两类以上“性质不同”的观测值赫尔默特方差分量估计就大概率比固定先验权更靠谱。2. 赫尔默特方差分量估计的核心原理2.1 残差二次型与方差分量的关系赫尔默特估计的基本思路并不复杂但推导过程容易劝退人我先尽量用直白的话把核心逻辑讲清楚。假设我们有 m 类观测值总体观测方程可以写成$$L A X \Delta$$其中 L 是 n 维观测向量A 是设计矩阵X 是待估参数向量Δ 是观测误差。将观测值按类别分块后每一类对应的设计矩阵为 A_i权阵为 P_i且假设各观测量之间不相关同一类内部可以相关赫尔默特基本形式按不相关处理相关观测值可以用LS-VCE方法扩展。对给定的一组权阵 P blockdiag(P_1, ..., P_m) 做加权最小二乘平差后可以得到残差向量 V。将残差按观测类别分块得到每类残差 V_i于是可以构造每类观测的残差二次型$$S_i V_i^T P_i V_i$$这个 S_i 乍看只是一个数却携带了关键信息它本质上是第 i 类观测噪声经平差投影后的“能量”。如果第 i 类观测的实际方差很大但权阵给得偏大即过于相信它那么平差会把这类观测的噪声强行压小残差反而会显得异常。反过来也一样。赫尔默特的核心数学结论是S_i 的期望值可以写成各类未知方差分量 σ_j² 的线性组合$$E(S_i) \sum_{j1}^{m} H_{ij} \sigma_j^2$$也就是说我们可以用观测到的 S 向量由 S_i 组成、通过系数矩阵 H反解出未知的方差分量 σ²。这个“反解”过程就是赫尔默特方差分量估计。2.2 估计矩阵H的构成与推导要点很多教程直接把 H 矩阵的公式甩出来却没有解释每个项的意义导致初学者即使会用代码也心里发虚。这里我把关键公式写出来并逐个拆解。先定义几个基础量。令整体法方程矩阵为$$N A^T P A \sum_{j1}^{m} A_j^T P_j A_j \sum_{j1}^{m} N_j$$其中 N_j 是第 j 类观测对法方程的贡献。再令 n_j 是第 j 类观测值的个数。赫尔默特估计矩阵 H 的每个元素为$$H_{ii} n_i - 2,\mathrm{tr}(N^{-1} N_i) \mathrm{tr}(N^{-1} N_i N^{-1} N_i)$$$$H_{ij} \mathrm{tr}(N^{-1} N_j N^{-1} N_i), \quad (i eq j)$$结合 S 向量方差分量的估计值就是$$\hat{\sigma}^2 H^{-1} S$$这几个项怎么理解n_i 是第 i 类观测的个数代表“完全自由”时的信息量。减去 2 tr(N^{-1}N_i) 是因为平差过程中引入了未知参数估计消耗了自由度每一个参数估计都会把一部分观测噪声“吸收”进参数里使得残差的能量小于原始观测噪声的能量。最后加上的 tr(N^{-1}N_iN^{-1}N_i) 是一个自相关修正项刻画第 i 类观测在得知自身噪声信息后对参数估计的二次反馈。非对角项 tr(N^{-1}N_jN^{-1}N_i) 则衡量第 j 类和第 i 类观测在参数估计中的相互纠缠程度——当两类观测都强约束同一个参数时这个项就会比较大。实际工程中我们一般不需要手动计算这些迹的值numpy或者MATLAB里一行 trace 就搞定了。但理解这些项的含义对排查异常结果非常有帮助。2.3 迭代流程与收敛逻辑赫尔默特方差分量估计不是一次计算就能完成的因为 S_i 是在给定权阵的前提下得到的残差二次型而权阵本身又要靠估计出的方差分量来修正。这个过程天然需要迭代。完整的迭代流程如下初始化给每类观测设定初始权阵 P_i。可以是等权也可以是高度角模型给出的先验权。工程上建议给一个相对合理但不追求精度的初值。做一次加权最小二乘平差得到参数估计值 X̂ 和残差 V并按观测类别切分得到 V_i。对每一类观测计算 S_i V_i^T P_i V_i。根据设计矩阵、当前权阵和法方程矩阵 N构造 H 矩阵。解线性方程 H σ² S得到方差分量估计值。更新权阵通常取第一类观测的方差分量为基准将各类观测的方差分量做归一化然后按 P_i ← P_i / σ̂_i² 更新权阵。重复步骤2到6直到归一化后的方差分量都接近1或者σ̂的变化量小于设定阈值即各类观测的实际噪声与当前权阵已经匹配。输出最后一次平差得到的参数估计值。收敛逻辑很直观当各类方差分量估计值都接近同一个尺度时说明当前权阵已经和实际噪声水平一致了继续迭代也不会带来明显改善可以收手。正常情况下3到10轮迭代就能收敛我在实际项目中很少见到超过15轮的。3. Python代码实现与模拟验证3.1 代码结构与关键函数说明代码我用Python numpy实现没有依赖任何GNSS专业库目的是让原理部分和代码一一对应方便读者对照理解。整个实现封装成一个类 HelmertVCE核心方法只有两个一个是组装整体权阵一个是执行单轮方差分量估计。我刻意把类写得简洁可读没有做太多防御性编程实际工程使用时可自行加参数校验和异常处理。代码的设计思路是输入按类别组织的设计矩阵列表、观测值列表、初始权阵列表然后调用 fit 方法完成迭代。3.2 模拟三类观测值的完整示例为了验证代码正确性我构造一个简单的线性回归问题3个未知参数3类观测值每类观测的真实噪声标准差不同。第一类标准差0.1第二类0.5第三类1.0。由于方差分量与标准差的平方成正比赫尔默特估计出来的相对方差分量理论上应接近 1:25:100。以下是完整代码可以直接复制运行import numpy as np class HelmertVCE: def __init__(self, A_list, L_list, P_list): self.A_list [] self.L_list [] self.P_list [] for A, L, P in zip(A_list, L_list, P_list): self.A_list.append(np.asarray(A, dtypefloat)) self.L_list.append(np.asarray(L, dtypefloat).reshape(-1, 1)) self.P_list.append(np.asarray(P, dtypefloat)) self.m len(A_list) # 观测类别数 self.n self.A_list[0].shape[1] # 待估参数个数 def _assemble(self): 组装整体设计矩阵、观测向量和分块对角权阵 A np.vstack(self.A_list) L np.vstack(self.L_list) P np.zeros((A.shape[0], A.shape[0])) start 0 for Pi in self.P_list: ni Pi.shape[0] P[start:start ni, start:start ni] Pi start ni return A, L, P def _vce_once(self, V_list, N): 单轮赫尔默特方差分量估计返回H矩阵和S向量 m self.m N_inv np.linalg.inv(N) H np.zeros((m, m)) S np.zeros(m) N_list [] for Ai, Pi in zip(self.A_list, self.P_list): N_list.append(Ai.T Pi Ai) for i in range(m): Vi V_list[i] Pi self.P_list[i] S[i] float(Vi.T Pi Vi) Ninv_Ni N_inv N_list[i] H[i, i] Pi.shape[0] - 2.0 * np.trace(Ninv_Ni) np.trace(Ninv_Ni Ninv_Ni) for k in range(m): if k i: continue H[i, k] np.trace(N_inv N_list[k] N_inv N_list[i]) return H, S def fit(self, max_iter50, tol1e-6): 迭代求解方差分量返回最终参数解、验后单位权方差和迭代历史 history [] for it in range(1, max_iter 1): A, L, P self._assemble() N A.T P A W A.T P L X np.linalg.solve(N, W) V A X - L V_list [] start 0 for Pi in self.P_list: ni Pi.shape[0] V_list.append(V[start:start ni]) start ni H, S self._vce_once(V_list, N) sigma2 np.linalg.solve(H, S) # 以第一类观测的方差分量为基准做归一化更新权阵 scale sigma2[0] if scale 0: scale 1.0 # 负方差处理正文4.1节会详细说明 sigma2_rel sigma2 / scale history.append(sigma2_rel.copy()) for j in range(self.m): self.P_list[j] self.P_list[j] / sigma2_rel[j] if np.max(np.abs(sigma2_rel - 1.0)) tol: break # 收敛后做最后一次平差返回最终结果 A, L, P self._assemble() N A.T P A W A.T P L X np.linalg.solve(N, W) V A X - L deg_free L.shape[0] - self.n sigma0_2 float(V.T P V / deg_free) return X, sigma0_2, history if __name__ __main__: rng np.random.default_rng(42) n_params 3 X_true np.array([3.0, -1.0, 2.0]) # 第一类20个观测噪声标准差0.1 A1 rng.normal(size(20, n_params)) L1 A1 X_true 0.1 * rng.normal(size20) P1 np.eye(20) # 第二类15个观测噪声标准差0.5 A2 rng.normal(size(15, n_params)) L2 A2 X_true 0.5 * rng.normal(size15) P2 np.eye(15) # 第三类10个观测噪声标准差1.0 A3 rng.normal(size(10, n_params)) L3 A3 X_true 1.0 * rng.normal(size10) P3 np.eye(10) # 赫尔默特方差分量估计 vce HelmertVCE([A1, A2, A3], [L1, L2, L3], [P1, P2, P3]) X_vce, sigma0, history vce.fit() print(真值: , X_true) print(VCE参数解: , X_vce.ravel()) print(估计的相对方差分量: , history[-1]) print(迭代轮数: , len(history)) # 对比直接用等权最小二乘不做VCE A_all np.vstack([A1, A2, A3]) L_all np.vstack([L1, L2, L3]) P_all np.eye(A_all.shape[0]) X_ols np.linalg.solve(A_all.T P_all A_all, A_all.T P_all L_all) print(等权参数解: , X_ols.ravel())3.3 运行结果解读与效果分析在这个模拟设置下理论上的相对方差分量是 1:25:100。运行代码后赫尔默特迭代得到的估计值会在理论值附近波动不会完全精确相等因为观测噪声是随机生成的有限样本下的估计本身带有随机性。但几个关键现象一定会出现第一迭代收敛很快通常5轮以内就能让归一化方差分量接近1。如果初始权阵给得离谱迭代轮数会略多但整体收敛趋势是稳定的。第二等权最小二乘得到的参数解会被第二类和第三类大噪声观测明显拉偏。第三类观测噪声标准差是第一类的10倍但等权平差却把它们当作同样可靠的信息源参数解自然偏离真值。赫尔默特估计在迭代中自动降低了第二类和第三类的权参数解会明显向真值靠拢。第三算出的验后单位权方差σ0²会比较接近1。这是因为收敛后各类观测的权与方差分量已经匹配单位权方差归一到合理范围。如果σ0²远大于1说明整体观测噪声被低估了权阵整体偏大如果远小于1则说明权阵整体偏小。我想强调一下代码里归一化这个细节。每次迭代解出的σ²是一个相对值它可以整体乘一个缩放因子而不影响参数解——因为参数估计只取决于各类观测权之间的比值。代码中选择固定第一类观测的方差为1相当于把第一类观测作为基准尺度其他类的权都相对于第一类进行调整。这个处理方式保证了迭代的数值稳定性。4. 工程应用中的常见问题与避坑指南4.1 方差分量为负怎么办赫尔默特估计最常见的“翻车现场”就是解出来的某个方差分量为负。方差按定义必须是非负的出现负值说明估计过程出了问题这时候不能简单地取绝对值糊弄过去要系统排查原因。第一个原因是先验权阵严重偏离实际。初始权如果给得太离谱第一轮迭代的残差二次型S就会携带很强的误导性信息H矩阵又基于当前权计算两者不匹配就可能解出负方差。解决方法是换一个更合理的初始权或者先跑几轮固定权重的预平差等残差稳定后再启动VCE。第二个原因是函数模型本身有误。比如有未发现的粗差、未估计的系统性偏差、未建模的对流层残差这些误差会全部“伪装”成观测噪声污染方差分量估计。这时即使迭代能跑通估计出的方差分量也不可信。我遇到过类似的情况某类观测值里混入了几个粗差方差分量估计结果直接崩溃删掉粗差后一切正常。第三个原因是某类观测值的数量太少或者这类观测对应的多余观测数不足。赫尔默特估计本质是从残差中提取信息如果某类观测只有两三个统计信息太少估计结果方差很大甚至出现负值。这时可以考虑把相近的观测类型合并成一个大类增加样本量。工程上的务实做法是出现负方差时先停下来检查函数模型和数据质量而不是简单地把负值替换成小正数继续迭代。4.2 迭代不收敛或振荡怎么办赫尔默特迭代大多数情况下很稳定但在某些特定条件下会出现不收敛或者来回振荡的情况。我自己总结下来最常见的原因有两个。一个是H矩阵病态。当两类观测高度相关、对参数的约束方式几乎相同时H矩阵中对应的行会接近线性相关导致 H 接近奇异解出的方差分量剧烈跳动。检查方法很简单看H矩阵的条件数如果 cond(H) 超过 1e12 基本就是病态了。缓解方法包括合并强相关的观测类、增加正则化项、或者换用最小二乘方差分量估计LS-VCE这类更稳健的方法。另一个是初始权阵太差导致迭代路径振荡。这时可以给权更新加一个阻尼因子例如每次迭代不是直接更新到目标值而是按 50% 的步长逼近$$P_i^{(k1)} P_i^{(k)} \left( \alpha / \sigma_i^2 (1-\alpha) / \bar{\sigma}^2 \right)$$这个式子看起来复杂实际效果就是让权阵“小步慢跑”牺牲一点收敛速度换取稳定性。我在处理多系统融合数据时用这个小技巧解决过不止一次振荡问题。另外提醒一点收敛阈值不要设得太苛刻。GNSS数据本身噪声大方差分量估计的随机波动是客观存在的tol 设为 1e-4 到 1e-6 就足够了再小没有实际意义只会增加无谓的迭代轮数。4.3 大型GNSS网中的计算效率问题上面给出的代码直接对法方程矩阵求逆这在参数个数几百上千的GNSS网平差中会非常吃力。实际工程中需要做几个方面的优化。第一避免显式求 N^{-1}。GNSS网平差的法方程矩阵通常是稀疏的尤其是站坐标、模糊度参数成百上千时直接求逆既慢又不稳定。正确做法是用Cholesky分解或者LDL^T分解把求逆转化为回代求解。第二迹的计算可以走捷径。H矩阵中需要计算 tr(N^{-1}N_i) 和 tr(N^{-1}N_iN^{-1}N_j)这两项如果直接按矩阵乘法算涉及大型矩阵乘法开销不小。对于分块对角结构的权阵可以利用“迹只取对角元素”的性质只计算 N^{-1}N_i 的对角元素避免完整矩阵相乘。第三观测类型划分不要过细。多系统融合时常规做法是把观测分为码伪距和载波相位两大类每个系统每一类再单独区分这样就有八个类别以上。从原理上说没问题但从统计稳健性角度讲类别越多每类样本量越少估计越不稳定。我的建议是优先按噪声水平差异最大的维度分比如伪距/相位、系统间高度角的影响通过先验权而不是VCE处理。第四迭代过程中的敏感度分析。大型网一旦出现异常每次完整迭代代价都不低建议在迭代初期就打印每轮方差分量的变化轨迹提前发现振荡或者负方差问题不要等到迭代完才发现白跑了一轮。4.4 赫尔默特与其他方差估计方法的关系说到赫尔默特很多人会提到MINQUE和LS-VCE这里顺带把关系理清楚避免大家混淆。MINQUE最小范数二次无偏估计Minimum Norm Quadratic Unbiased Estimation是Rao提出的一套系统化的方差分量估计理论。它的显著特点是不需要迭代一次性给出方差分量的估计值但需要预先设定先验方差比。赫尔默特可以看作是MINQUE的迭代化实现每一步都用当前估计出的方差分量更新权阵再重新估计直到收敛。LS-VCELeast-Squares Variance Component Estimation是Teunissen提出的一种更现代的框架它的核心思想是把方差分量估计问题本身也构造成一个最小二乘问题。相比赫尔默特LS-VCE更容易处理观测值之间存在相关性的情况而且可以很方便地引入方差分量的先验信息。代价是数学上更复杂计算量也更大。在实际GNSS数据处理项目中我的排序是观测值不相关的经典场景直接上赫尔默特简单高效观测值相关性强或者对统计性质要求高的研究场景优先考虑LS-VCEMINQUE更多作为一种理论工具在软件实现中其实用得不多。最后分享一个我自己的使用习惯。在项目里我不会一上来就开赫尔默特而是先用高度角模型给一个不差的初值再跑两三轮VCE这样既快又稳直接给等权虽然也能收敛但前期迭代容易抖动。另外我会把每一轮各类观测的权更新日志打出来一旦发现某两类权在反复横跳基本可以断定函数模型或者粗差处理出了问题——这时候修VCE没用回头查数据更实在。赫尔默特不是万能的但它确实是让随机模型从“拍脑袋”走向“用数据说话”的可靠工具。