有一次外业收工回来六台接收机、三个同步时段基线解算软件吐出上百条 ΔX、ΔY、ΔZ每条基线后面还挂着一个 3×3 的方差-协方差阵。我当时的想法很朴素把基线读进来点一下平差坐标就出来了。结果第一次运行直接卡住软件用红字圈出一个闭合环闭合差 11 ppm不让我往下走。后来才慢慢明白GPS网平差真正难的地方根本不在最小二乘那几行公式而在前面挂着的一整套质量检验、中间夹着的基准与随机模型取舍、后面接着的坐标转换和高程处理。任何一环理解偏了解出来的坐标都会看起来特别漂亮——单位权中误差 1.000、改正数全是零点几毫米但跟真实位置差出十几厘米而且你还找不出错在哪。这篇内容面向的是手上有真实观测数据、需要自己把控平差过程的人做控制网的技术员、写 GNSS 数据处理脚本的开发者、以及那些用现成软件但总觉得哪里不对的使用者。我会把闭合环检验、间接平差的函数模型与随机模型、基准定义、Python 实现、约束平差与坐标系转换、大地高与地图偏移这几件事串成一条线讲清楚中间穿插我自己踩过、也确实让我返过工的细节。1. 闭合环先跑一遍平差软件为什么宁可罢工也不出结果1.1 基线向量是差不是坐标相对定位解出来的东西本质上就是两台接收机之间的三维坐标差ΔX_ij X_j − X_i ε它天生没有绝对位置。一条基线只告诉你从 A 到 B 走了多远、朝哪个方向至于 A 本身在哪基线本身一个字都没说。要得到绝对坐标必须人为指定一个基准——固定某个点的坐标或者固定网的重心然后靠这些差把其余点推出来。很多人第一次做平差时会困惑为什么未知数个数是点数乘 3而观测方程个数是基线数乘 3两者对不上还要求解原因就在这——没有基准方程组的系数矩阵天然秩亏 3解出来的是一族无穷多的解只是形状相同。理解了这一点闭合环检验的意义就清楚了。三条基线 A→B、B→C、C→A 围成一个三角形理论上绕一圈必须精确回到原点也就是V_ΔX ΔX_AB ΔX_BC ΔX_CA ≈ 0实际观测永远不严格为零这个残差就是环闭合差。它是完全独立于基准的内部一致性度量——换什么基准、怎么约束闭合差都不会变。所以行业里把它放在平差之前当作准入门槛环都闭不上后面的平差只是在给错误数据做精致包装。1.2 同步环、异步环、重复基线三道不同的筛子这三类检验经常被混为一谈但它们的严格程度和指示的问题完全不同。同步环是同一观测时段内多台接收机同时采集、由基线连成的闭合环因为这些基线共享同一批卫星、同一个时段的大气条件误差高度相关闭合差通常很小所以它的限值最严。异步环由不同时段的基线拼成误差来源独立叠加闭合差会明显大一些。重复基线则是对同一条边的多次观测结果做比较直接反映这条边自身的稳定性。检验类型观测构成主要暴露的问题常见限值参考同步环同一时段内多条基线闭合解算软件逻辑错误、点名映射错误分量闭合差 ≤ √(3n)·σ全长相对闭合差按等级可达 1~2 ppm异步环跨时段基线闭合大气、星历、多路径、粗差分量闭合差 ≤ 3√n·σ相对闭合差通常 2~5 ppm重复基线同一条边多次解算单条基线的可靠性较差 ≤ 2√2·σ注意这里的 σ 指的是基线解算给出的该边中误差n 是环上的边数。不同行业规范控制测量、变形监测、工程测量给的系数并不一致动手前先翻你手上那份规范别照搬网上的数字。1.3 闭合差超限时的排查顺序超限之后最常见的反应是去调解算参数、改截止高度角、换星历试图把闭合差压下去。我的经验是先按下面这个顺序走八成问题在第一、二步就能找到第一检查点名映射。这是最容易被忽略又最致命的一类问题。多台接收机、多时段观测时如果某个文件里的点名重复或者仪器内点位编号和记录簿对不上软件读进来的基线就变成了另一张网。表现往往就是某几段闭合差诡异超限其他环却干干净净。第二核对天线高。天线高量错 1 cm在 500 m 短边上就是 20 ppm在 5 km 边上才 2 ppm。量高错误典型特征是短边闭合差特别大、长边正常而且超限的环全部与同一个点有关。这个规律记住能省你半天时间。第三看单条基线的解算质量。RMS、Ratio、数据利用率、模糊度固定情况哪条边明显差先把它剔掉重算看闭合差有没有改善。第四看卫星几何与观测时段。PDOP 长时间的时段、卫星数不足的时段、多路径严重的站点比如旁边有大面积水面或玻璃幕墙解出来的基线系统性偏差会比较大。第五考虑星历。长基线几十公里以上用广播星历和用精密星历差异可以到 ppm 量级短基线几乎没区别。基线一长这件事就不能省。2. 间接平差的骨架每一项到底对应现实里的什么东西2.1 观测方程 V Bδx̂ − l 的逐项拆解网平差绝大多数情况下用间接平差也叫参数平差。写成矩阵形式就是V B·δx̂ − l这四个符号对应到 GPS 网里是这样l常数项向量观测值减去用近似坐标算出来的观测值。对一条基线 i→j 的 X 分量l ΔX_ij^观测 − (X_j⁰ − X_i⁰)其中 X⁰ 是近似坐标。近似坐标的精度要求不高跟真实坐标差个几百米都不影响最终解因为一阶偏导是 ±1线性化误差可以忽略。这也是为什么 GPS 网平差不需要迭代。δx̂未知数向量每个待定点的三维坐标改正数按点号顺序排成 3n×1 的列向量。B设计矩阵/系数矩阵每一行对应一个观测分量每一列对应一个未知数。对基线 i→j 的 X 分量那一行只有 i 点的 X 列是 −1j 点的 X 列是 1其余全是 0。整张 B 矩阵是极其稀疏的稀疏度一般在 99% 以上。V改正数向量观测值的平差改正数是拿来判断有没有粗差的关键指标。法方程就是 N·δx̂ W其中 N BᵀPBW BᵀPl。看到这里你会发现整件事的核心工作量不在于解方程而在于怎么构造 P权阵和怎么处理 N 的秩亏。2.2 随机模型那份几乎没人认真读过的协方差阵基线解算软件输出的不只是 ΔX、ΔY、ΔZ还有每条基线对应的 3×3 协方差阵。里面包含了三个分量各自的方差以及它们两两之间的协方差。这一点很重要同一基线里 ΔX、ΔY、ΔZ 的误差不是独立的因为它们在同一个解算过程里被一起估计出来尤其在长基线上几十公里相关性很强。权阵的构造方式是 P σ₀²·D⁻¹。D 是全部基线的协方差阵如果只取独立基线D 可以写成按基线分块的对角块结构求逆非常快如果用了全部基线基线之间强相关D 就不再是块对角这个坑在下一节展开。实践中最常见的做法是引入一个比例因子scale factor把所有基线的协方差阵统一乘以同一个 k然后调 k让平差出来的单位权中误差接近 1。这个做法的合理性在于——基线解算给出的协方差阵通常偏乐观而且没有反映网内部的真实符合程度。k 的取值通常在 1 到 100 之间。我的习惯是先用 k1 跑一次自由网看单位权中误差的量级再反推一个 k 的初值。一个非常隐蔽的坑不同软件输出协方差阵的单位不一样。有的写 mm²有的写 m²有的把标准差当方差写进去。单位差 10⁶ 倍的时候虽然最终的坐标解不受影响整体缩放但单位权中误差和点位精度会全错卡方检验也会失效。动手写解析脚本前务必先拿一条已知边手工验算一下量级。2.3 基准怎么定自由网、固定基准、拟稳平差的取舍基准问题决定了解相对于什么。三种主流选择自由网最小约束平差只施加恰好消除秩亏的最少约束——通常是固定网的重心或者用伪逆求最小范数解。它得到的是网的纯形状专门用来检核内部符合精度、发现粗差、以及与已知点做兼容性对比。它不给你绝对坐标但它是所有后续步骤的基础。固定基准平差固定一个点的三维坐标再加上一个方位约束。这种方案简单但把全部基准压力都压在一个点上——这个点但凡有点位移或坐标误差整个网跟着转。拟稳平差先把网中相对稳定的点挑出来判断依据通常是自由网解里的位移量以这些点为基准。变形监测网基本都用它工程控制网在已知点质量参差时也值得用。约束平差则是固定多个已知点把 GPS 网的基准强行对齐到地方坐标系。它的问题在于已知点本身有误差约束进去就是硬把误差焊死。所以规范流程一般是三维无约束平差 → 三维约束平差顺便做已知点兼容性检验→ 二维约束平差得到地方平面坐标。中间那步兼容性检验千万别跳。3. 用 Python 把一次三维无约束平差完整跑通3.1 基线文件的组织方式与读取我用得最多的中间格式是自己定义的一份纯文本每行一条基线后面跟六个协方差参数三个方差 三个协方差。这样无论是从软件导出的报告里解析还是从解算程序直接输出格式都很干净# from to dX dY dZ var_XX var_YY var_ZZ cov_XY cov_XZ cov_YZ A B 1234.5678 -2345.6789 3456.7890 4.0 4.5 5.2 0.3 0.2 0.4 A C -3456.1234 1234.5678 -2345.6789 5.1 5.5 6.0 0.4 0.3 0.5 ...单位统一用米和平方米。读取的时候用 numpy顺手把点名映射成索引import numpy as np def load_baselines(path): rows [] with open(path) as f: for line in f: line line.strip() if not line or line.startswith(#): continue parts line.split() rows.append(parts) return rows def build_index(rows): pts [] for r in rows: for p in (r[0], r[1]): if p not in pts: pts.append(p) return pts, {p: i for i, p in enumerate(pts)}3.2 组 B矩阵符号和索引是最容易翻车的地方B 矩阵的构造逻辑简单但写错一次就会得到完全不合理的结果而且错误不会报异常——它会给你一个看起来正常的解。所以每次写完这段代码我都会拿一条基线手工验算一遍。def build_B(rows, idx, n_pts): n_obs 3 * len(rows) n_unk 3 * n_pts B np.zeros((n_obs, n_unk)) for k, r in enumerate(rows): i idx[r[0]] j idx[r[1]] for c in range(3): B[3*k c, 3*i c] -1.0 # 起点负 B[3*k c, 3*j c] 1.0 # 终点正 return B符号约定必须和常数项 l 的计算方式严格一致。如果 l 定义成观测减近似那么起点就是 −1、终点 1如果反过来定义成近似减观测符号要全部对调否则解出来的坐标改正数方向是反的最终坐标会错得离谱。另一个常见错误是索引错位B 矩阵的第 3k、3k1、3k2 行分别对应 X、Y、Z如果列方向上也用同样的排列3i、3i1、3i2一一对应即可但一旦有人在某处把点号和分量的顺序写混结果就是 X 分量的观测被算到了 Y 的未知数上。3.3 秩亏与基准约束的两种处理方式自由网的 N 矩阵秩亏 3直接求逆会报错或者给出垃圾解。两种处理方式我都用过第一种是用最小二乘的通用求解器拿最小范数解dx, *_ np.linalg.lstsq(B, l, rcondNone)它的数学含义是求满足观测方程且改正数平方和最小的那个解同时范数最小等价于固定网的重心。做内部检核很方便缺点是基准的物理含义不直观。第二种是显式施加约束构造扩展法方程# C 是约束矩阵例如固定 0 号点的三个坐标 C np.zeros((3, 3 * n_pts)) C[0, 0] C[1, 1] C[2, 2] 1.0 N B.T P B W B.T P l A np.block([[N, C.T], [C, np.zeros((3, 3))]]) sol np.linalg.solve(A, np.concatenate([W, np.zeros(3)])) dx sol[:3 * n_pts]这种做法的好处是约束的物理意义明确我确实要固定这个点而且精度评定可以直接从 N 的逆里取。坏处是如果 C 的构造有问题法方程可能依然奇异numpy 会抛 LinAlgError这时候要回头检查 C 的秩是不是 3。3.4 精度评定单位权中误差、点位误差、卡方检验解完之后不做精度评定等于白做。三个必看的量V B dx - l sigma0 np.sqrt(V.T P V / (3 * len(rows) - (3 * n_pts - 3))) Qxx np.linalg.inv(N) point_std [sigma0 * np.sqrt(Qxx[3*i, 3*i] Qxx[3*i1, 3*i1] Qxx[3*i2, 3*i2]) for i in range(n_pts)]单位权中误差 σ₀是整体拟合质量的指标。它接近 1说明协方差阵的量级和实际符合程度匹配远大于 1说明观测精度没有协方差阵声称的那么好或者附近有粗差远小于 1说明协方差阵偏保守或者人为加了比例因子。它大于 1 是常态不必惊慌。点位中误差由 Qxx 的对角块给出。注意 Qxx 是加了约束之后的逆自由网里讨论点位精度本身意义不大因为那取决于你怎么定基准。卡方检验用来判断模型是否成立VᵀPV / σ₀² 应服从自由度等于多余观测数的 χ² 分布。落在两侧的临界值之外就要考虑随机模型是否设置错了。三个指标里我最看重的是改正数的分布。把所有基线的 V 按分量画出来看如果只是随机散布在零附近、没有超过 3σ 的心里就稳了如果出现某几条边的改正数同号且偏大通常是粗差或者这些边和其他边不是同一批数据解出来的比如混了不同星历、不同时段的高度角设置。4. 从 WGS-84 落到地方坐标约束平差与转换参数的顺序问题4.1 三维无约束到二维约束中间不能跳步完整的规范流程是三步走我见过太多人直接从三维基线跳到平面坐标结果一堆莫名其妙的问题第一步三维无约束平差。目的是检核网本身的内部一致性、剔除粗差、评估相对精度。这一步的成果是相对于网重心的坐标不能直接拿去用。第二步三维约束平差。选部分已知点参与约束把网的基准对齐到目标坐标系比如 CGCS2000 或地方坐标系的空间直角坐标。这一步的关键产出是已知点兼容性检验——比较已知点坐标和约束后坐标的残差看哪些点不兼容。不兼容的点要么剔除要么降权。第三步二维约束平差。在平面坐标系里做得到最终的地方平面坐标。这一步会引入投影和方位必须先把已知点的高斯平面坐标反算到椭球面或者直接在三维空间里做平面约束。跳过第二步的风险在于你把所有已知点一股脑约束进去其中某个点有 5 cm 的误差整个网会被它拖着转而你在平面坐标里看不出任何异常。4.2 七参数与四参数选错就差一个数量级布尔沙七参数模型三个平移、三个旋转、一个尺度共七个参数。适合大范围、跨投影带、三维坐标之间的转换比如 WGS-84 到 CGCS2000 或到地方独立坐标系。平面四参数两个平移、一个旋转、一个尺度共四个参数。适合小范围、同一投影带内、二维平面坐标之间的转换比如两套城市坐标之间的对接。选择依据就两条控制点的空间分布范围和坐标系基准是否一致。控制点范围超过一个投影带大约 3° 经差或者跨省用七参数四参数在大范围里会把投影变形吃进旋转参数里越远偏得越狠。控制点集中在一个城市内、都在同一带四参数足够而且四参数需要的最少公共点更少。高差大、点位分布跨越较大高程带必须用七参数因为四参数在平面上没法处理高程方向的影响。参数的个数直接决定了最少公共点数七参数至少 3 个四参数至少 2 个。但这只是数学上的最低要求实际使用中公共点数量应该是参数个数的两倍以上否则参数解没有冗余、无法检验一个点有粗差就能把整套参数带偏。解出参数后一定要留出几个点做检核看残差是否在厘米级。4.3 边长反算与投影改化这两项不做网就永远对不上这是最容易被跳过、后果也最严重的一步。GPS 基线是从空间得到的斜距或椭球面距离要拿去和地面网比较或者输入平面坐标必须做两项归算第一项高程归算归算到参考椭球面。近似公式是 ΔS₁ ≈ −Δh²/(2S)其中 Δh 是相对于参考椭球面的高差。举个数平均高程 100 m边长 1 km改正量约 −5 mm/km。听着不大但如果测区平均高程 500 m改正量就变成 −125 mm/km也就是 1 km 差 12.5 cm这已经远超多数工程网的允许误差了。第二项高斯投影改化归算到高斯平面。近似公式是 ΔS₂ ≈ S·y²/(2R²)y 是测区到中央子午线的平均横坐标R 是地球平均曲率半径。在 y 50 km 处投影变形约 1/32000也就是 1 km 差 3 cmy 100 km 处1 km 差 13 cm。两项的符号通常相反会有一部分抵消但绝不能指望它们自己抵完。判断测区是否要做投影改化的经验值是相对变形的绝对值超过 1/40000就必须处理。碰到投影变形过大的测区常见的工程解法是建立地方独立坐标系——选一个合适的抵偿面或者中央子午线把变形控制在可接受范围内。一个真实的返工案例测区平均高程 800 m、距中央子午线 60 km两项改正叠加后边长相对误差约 1/25000。当时没做归算直接拿 GPS 坐标算出来的边和地面网差了几十厘米第一反应以为是 GPS 数据有问题查了两天闭合环、换了三次星历最后发现是根本没做投影改化。5. 大地高、正常高以及把点画到地图上的那些偏移5.1 大地高不是海拔差的是高程异常GPS 直接解算出来的是大地高 H也就是点到参考椭球面的距离。我们平时说的海拔是正常高 h是点到似大地水准面的距离。两者之间的关系是H h ζζ 就是高程异常。它的量级在全国范围内从几米到几十米不等跟地形和地壳结构有关。如果你直接把 GPS 解出的大地高当成海拔填进成果表误差就是几十米这个量级——这个错误不致命因为太好发现了但它会让后续所有基于高程的处理都错位。5.2 GPS 高程拟合能用到什么范围把 GPS 和水准结合用少量重合点拟合出测区的高程异常曲面从而把大量 GPS 点的大地高换算成正常高这是常规做法。有用的经验是二次曲面拟合至少需要 6 个均匀分布的重合点。低于这个数曲面系数解不出来或者极不稳定。拟合范围要覆盖所有待求点不要外推。外推超过重合点包络范围 10% 以外精度下降得非常快。地形起伏大的地方山地区曲面拟合基本没法用。这时改用高程异常模型比如 EGM2008 之类的重力场模型给出一个初值精度大概在分米级对一般工程够用。平原地区、重合点分布均匀、点距 2~5 km二次曲面拟合的精度通常能到 2~3 cm这是比较乐观的情况。还有一个细节水准点的高程本身可能有沉降或者历史遗留的误差。拟合前先做一次重合点之间的高程异常变化趋势检查看有没有明显的异常点别把脏数据喂进拟合。5.3 BLH 到高德坐标GCJ-02 不是一个平移量数据要落到地图上就会碰到坐标系的问题。目前常见的有三套坐标系使用场景与 WGS-84 的关系WGS-84GPS 原始输出、国际标准基准GCJ-02高德、腾讯等国内地图非线性偏移偏移量随位置变化BD-09百度地图在 GCJ-02 基础上再做一次偏移关键点是 GCJ-02 的偏移不是简单加一个常数而是一组基于经纬度的非线性变换。所以我用某个点算出来的偏移量直接加在所有点上这种做法在几公里范围内大概能用跨几十公里就完全不准了。公开实现的思路大致是这样import math a 6378245.0 ee 0.00669342162296594323 def _tf_lat(lng, lat): ret -100.0 2.0*lng 3.0*lat 0.2*lat*lat 0.1*lng*lat 0.2*math.sqrt(abs(lng)) ret (20.0*math.sin(6.0*lng*math.pi) 20.0*math.sin(2.0*lng*math.pi)) * 2.0/3.0 ret (20.0*math.sin(lat*math.pi) 40.0*math.sin(lat/3.0*math.pi)) * 2.0/3.0 ret (160.0*math.sin(lat/12.0*math.pi) 320.0*math.sin(lat*math.pi/30.0)) * 2.0/3.0 return ret def _tf_lng(lng, lat): ret 300.0 lng 2.0*lat 0.1*lng*lng 0.1*lng*lat 0.1*math.sqrt(abs(lng)) ret (20.0*math.sin(6.0*lng*math.pi) 20.0*math.sin(2.0*lng*math.pi)) * 2.0/3.0 ret (20.0*math.sin(lng*math.pi) 40.0*math.sin(lng/3.0*math.pi)) * 2.0/3.0 ret (150.0*math.sin(lng/12.0*math.pi) 300.0*math.sin(lng/30.0*math.pi)) * 2.0/3.0 return ret def wgs84_to_gcj02(lng, lat): dlat _tf_lat(lng - 105.0, lat - 35.0) dlng _tf_lng(lng - 105.0, lat - 35.0) radlat lat / 180.0 * math.pi magic math.sin(radlat) magic 1 - ee * magic * magic sqrtmagic math.sqrt(magic) dlat (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi) dlng (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi) return lng dlng, lat dlat反方向GCJ-02 回 WGS-84没有解析解通常用迭代逼近或者一次近似。这部分有两个坑值得说一是高德的接口参数顺序是经度,纬度跟平时说经纬度的顺序是反的写错一次你会看到点跑到了另一个省二是这个转换本身有 1~2 m 量级的误差需要亚米级精度的成果绝对不能靠它得用带明确参数的坐标转换。成果导出到地图上做可视化时我习惯先转成 GeoJSON 或者带表头的 CSV把点号、转换后的经纬度、点位中误差一起输出。点位中误差直接作为属性字段带上在地图上按大小渲染颜色一眼就能看出网的精度分布——比看表格快得多。6. 复盘几个真踩过的坑6.1 单位和比例因子mm²、m²、微弧度单位问题我在前面提过一次但值得再强调因为它出现的频率实在太高。三种具体表现协方差阵用 mm² 而坐标用 m角度相关的量用微弧度而代码里按弧度算比例 ppm 和纯数字混用。这三种错误的共同特征是结果看起来合理但精度指标荒谬——坐标本身对不上是运气对上了是巧合。我现在的做法是写一个单位检查函数在读入数据后立刻验证某条已知边的中误差算出来应该在 5 mm 到 2 cm 之间如果算出来是 0.000005 或者 5000就说明单位有问题。这个检查花不了十分钟能省掉一整天。6.2 已知点兼容性检验不做等于把误差焊死约束平差里最危险的操作是把所有已知点都约束上去。已知点数量多的时候看起来约束越强、结果越可靠实际上是把每个已知点的误差都硬塞进了网形里。经典的判断指标是无约束平差坐标和已知坐标之间的差值如果某个已知点的差值明显大于其他点、且方向与周围点不一致它大概率有问题。处理的顺序是先用全部已知点做一次约束平差看残差分布把残差最大的点剔除重做重复直到剩下的点残差都在限差内。剔除的点不是数据坏了可能只是它属于另一套基准、或者经历了沉降但如果强行用最终成果的责任在你。一个具体的限差参考已知点残差超过 3 倍的点位中误差或者超过规范规定的已知点限差就应该被怀疑。不同等级控制网的具体数值差别很大别照搬。6.3 观测数据本身的可信度别指望平差兜底这是我最想强调的一点也是很多人忽略的平差是数据处理的最后一道关口它只能处理符合模型假设的误差。如果观测数据本身出了问题——信号受到干扰、载噪比在某个时段异常下降、伪距残差呈现整体偏移、多路径在特定方位角上持续存在——平差出来的坐标依然会是一个数学上完美的解单位权中误差可能还特别漂亮但它和真实位置的关系已经断了。所以在平差之前我会做几件事看每个站每个时段的载噪比时间序列异常下降或者分布异常的都拎出来看伪距和相位残差的空间分布如果残差在某个方位角上系统偏大多半是多路径那个方向的观测权重应该降低把同一时段的不同解算策略跑一遍做交叉验证结果差异大就说明数据本身不稳定。另外平差结果还可以用来反查数据问题如果改正数的空间分布呈现明显的系统性方向偏移比如所有基线在某个方向上改正数都是正的那多半不是随机误差而是数据处理链路上游的某个环节出了系统性问题。这种时候回头看原始观测比在平差环节调参数有用得多。说到底GPS 网平差教会我的不是那几行最小二乘公式而是先怀疑数据再怀疑模型最后才怀疑算法这个顺序。我见过太多人在算法里绕圈子其实问题一开始就在天线高、点名或者单位上。
GPS网平差全流程：闭合环检验、间接平差与坐标转换
有一次外业收工回来六台接收机、三个同步时段基线解算软件吐出上百条 ΔX、ΔY、ΔZ每条基线后面还挂着一个 3×3 的方差-协方差阵。我当时的想法很朴素把基线读进来点一下平差坐标就出来了。结果第一次运行直接卡住软件用红字圈出一个闭合环闭合差 11 ppm不让我往下走。后来才慢慢明白GPS网平差真正难的地方根本不在最小二乘那几行公式而在前面挂着的一整套质量检验、中间夹着的基准与随机模型取舍、后面接着的坐标转换和高程处理。任何一环理解偏了解出来的坐标都会看起来特别漂亮——单位权中误差 1.000、改正数全是零点几毫米但跟真实位置差出十几厘米而且你还找不出错在哪。这篇内容面向的是手上有真实观测数据、需要自己把控平差过程的人做控制网的技术员、写 GNSS 数据处理脚本的开发者、以及那些用现成软件但总觉得哪里不对的使用者。我会把闭合环检验、间接平差的函数模型与随机模型、基准定义、Python 实现、约束平差与坐标系转换、大地高与地图偏移这几件事串成一条线讲清楚中间穿插我自己踩过、也确实让我返过工的细节。1. 闭合环先跑一遍平差软件为什么宁可罢工也不出结果1.1 基线向量是差不是坐标相对定位解出来的东西本质上就是两台接收机之间的三维坐标差ΔX_ij X_j − X_i ε它天生没有绝对位置。一条基线只告诉你从 A 到 B 走了多远、朝哪个方向至于 A 本身在哪基线本身一个字都没说。要得到绝对坐标必须人为指定一个基准——固定某个点的坐标或者固定网的重心然后靠这些差把其余点推出来。很多人第一次做平差时会困惑为什么未知数个数是点数乘 3而观测方程个数是基线数乘 3两者对不上还要求解原因就在这——没有基准方程组的系数矩阵天然秩亏 3解出来的是一族无穷多的解只是形状相同。理解了这一点闭合环检验的意义就清楚了。三条基线 A→B、B→C、C→A 围成一个三角形理论上绕一圈必须精确回到原点也就是V_ΔX ΔX_AB ΔX_BC ΔX_CA ≈ 0实际观测永远不严格为零这个残差就是环闭合差。它是完全独立于基准的内部一致性度量——换什么基准、怎么约束闭合差都不会变。所以行业里把它放在平差之前当作准入门槛环都闭不上后面的平差只是在给错误数据做精致包装。1.2 同步环、异步环、重复基线三道不同的筛子这三类检验经常被混为一谈但它们的严格程度和指示的问题完全不同。同步环是同一观测时段内多台接收机同时采集、由基线连成的闭合环因为这些基线共享同一批卫星、同一个时段的大气条件误差高度相关闭合差通常很小所以它的限值最严。异步环由不同时段的基线拼成误差来源独立叠加闭合差会明显大一些。重复基线则是对同一条边的多次观测结果做比较直接反映这条边自身的稳定性。检验类型观测构成主要暴露的问题常见限值参考同步环同一时段内多条基线闭合解算软件逻辑错误、点名映射错误分量闭合差 ≤ √(3n)·σ全长相对闭合差按等级可达 1~2 ppm异步环跨时段基线闭合大气、星历、多路径、粗差分量闭合差 ≤ 3√n·σ相对闭合差通常 2~5 ppm重复基线同一条边多次解算单条基线的可靠性较差 ≤ 2√2·σ注意这里的 σ 指的是基线解算给出的该边中误差n 是环上的边数。不同行业规范控制测量、变形监测、工程测量给的系数并不一致动手前先翻你手上那份规范别照搬网上的数字。1.3 闭合差超限时的排查顺序超限之后最常见的反应是去调解算参数、改截止高度角、换星历试图把闭合差压下去。我的经验是先按下面这个顺序走八成问题在第一、二步就能找到第一检查点名映射。这是最容易被忽略又最致命的一类问题。多台接收机、多时段观测时如果某个文件里的点名重复或者仪器内点位编号和记录簿对不上软件读进来的基线就变成了另一张网。表现往往就是某几段闭合差诡异超限其他环却干干净净。第二核对天线高。天线高量错 1 cm在 500 m 短边上就是 20 ppm在 5 km 边上才 2 ppm。量高错误典型特征是短边闭合差特别大、长边正常而且超限的环全部与同一个点有关。这个规律记住能省你半天时间。第三看单条基线的解算质量。RMS、Ratio、数据利用率、模糊度固定情况哪条边明显差先把它剔掉重算看闭合差有没有改善。第四看卫星几何与观测时段。PDOP 长时间的时段、卫星数不足的时段、多路径严重的站点比如旁边有大面积水面或玻璃幕墙解出来的基线系统性偏差会比较大。第五考虑星历。长基线几十公里以上用广播星历和用精密星历差异可以到 ppm 量级短基线几乎没区别。基线一长这件事就不能省。2. 间接平差的骨架每一项到底对应现实里的什么东西2.1 观测方程 V Bδx̂ − l 的逐项拆解网平差绝大多数情况下用间接平差也叫参数平差。写成矩阵形式就是V B·δx̂ − l这四个符号对应到 GPS 网里是这样l常数项向量观测值减去用近似坐标算出来的观测值。对一条基线 i→j 的 X 分量l ΔX_ij^观测 − (X_j⁰ − X_i⁰)其中 X⁰ 是近似坐标。近似坐标的精度要求不高跟真实坐标差个几百米都不影响最终解因为一阶偏导是 ±1线性化误差可以忽略。这也是为什么 GPS 网平差不需要迭代。δx̂未知数向量每个待定点的三维坐标改正数按点号顺序排成 3n×1 的列向量。B设计矩阵/系数矩阵每一行对应一个观测分量每一列对应一个未知数。对基线 i→j 的 X 分量那一行只有 i 点的 X 列是 −1j 点的 X 列是 1其余全是 0。整张 B 矩阵是极其稀疏的稀疏度一般在 99% 以上。V改正数向量观测值的平差改正数是拿来判断有没有粗差的关键指标。法方程就是 N·δx̂ W其中 N BᵀPBW BᵀPl。看到这里你会发现整件事的核心工作量不在于解方程而在于怎么构造 P权阵和怎么处理 N 的秩亏。2.2 随机模型那份几乎没人认真读过的协方差阵基线解算软件输出的不只是 ΔX、ΔY、ΔZ还有每条基线对应的 3×3 协方差阵。里面包含了三个分量各自的方差以及它们两两之间的协方差。这一点很重要同一基线里 ΔX、ΔY、ΔZ 的误差不是独立的因为它们在同一个解算过程里被一起估计出来尤其在长基线上几十公里相关性很强。权阵的构造方式是 P σ₀²·D⁻¹。D 是全部基线的协方差阵如果只取独立基线D 可以写成按基线分块的对角块结构求逆非常快如果用了全部基线基线之间强相关D 就不再是块对角这个坑在下一节展开。实践中最常见的做法是引入一个比例因子scale factor把所有基线的协方差阵统一乘以同一个 k然后调 k让平差出来的单位权中误差接近 1。这个做法的合理性在于——基线解算给出的协方差阵通常偏乐观而且没有反映网内部的真实符合程度。k 的取值通常在 1 到 100 之间。我的习惯是先用 k1 跑一次自由网看单位权中误差的量级再反推一个 k 的初值。一个非常隐蔽的坑不同软件输出协方差阵的单位不一样。有的写 mm²有的写 m²有的把标准差当方差写进去。单位差 10⁶ 倍的时候虽然最终的坐标解不受影响整体缩放但单位权中误差和点位精度会全错卡方检验也会失效。动手写解析脚本前务必先拿一条已知边手工验算一下量级。2.3 基准怎么定自由网、固定基准、拟稳平差的取舍基准问题决定了解相对于什么。三种主流选择自由网最小约束平差只施加恰好消除秩亏的最少约束——通常是固定网的重心或者用伪逆求最小范数解。它得到的是网的纯形状专门用来检核内部符合精度、发现粗差、以及与已知点做兼容性对比。它不给你绝对坐标但它是所有后续步骤的基础。固定基准平差固定一个点的三维坐标再加上一个方位约束。这种方案简单但把全部基准压力都压在一个点上——这个点但凡有点位移或坐标误差整个网跟着转。拟稳平差先把网中相对稳定的点挑出来判断依据通常是自由网解里的位移量以这些点为基准。变形监测网基本都用它工程控制网在已知点质量参差时也值得用。约束平差则是固定多个已知点把 GPS 网的基准强行对齐到地方坐标系。它的问题在于已知点本身有误差约束进去就是硬把误差焊死。所以规范流程一般是三维无约束平差 → 三维约束平差顺便做已知点兼容性检验→ 二维约束平差得到地方平面坐标。中间那步兼容性检验千万别跳。3. 用 Python 把一次三维无约束平差完整跑通3.1 基线文件的组织方式与读取我用得最多的中间格式是自己定义的一份纯文本每行一条基线后面跟六个协方差参数三个方差 三个协方差。这样无论是从软件导出的报告里解析还是从解算程序直接输出格式都很干净# from to dX dY dZ var_XX var_YY var_ZZ cov_XY cov_XZ cov_YZ A B 1234.5678 -2345.6789 3456.7890 4.0 4.5 5.2 0.3 0.2 0.4 A C -3456.1234 1234.5678 -2345.6789 5.1 5.5 6.0 0.4 0.3 0.5 ...单位统一用米和平方米。读取的时候用 numpy顺手把点名映射成索引import numpy as np def load_baselines(path): rows [] with open(path) as f: for line in f: line line.strip() if not line or line.startswith(#): continue parts line.split() rows.append(parts) return rows def build_index(rows): pts [] for r in rows: for p in (r[0], r[1]): if p not in pts: pts.append(p) return pts, {p: i for i, p in enumerate(pts)}3.2 组 B矩阵符号和索引是最容易翻车的地方B 矩阵的构造逻辑简单但写错一次就会得到完全不合理的结果而且错误不会报异常——它会给你一个看起来正常的解。所以每次写完这段代码我都会拿一条基线手工验算一遍。def build_B(rows, idx, n_pts): n_obs 3 * len(rows) n_unk 3 * n_pts B np.zeros((n_obs, n_unk)) for k, r in enumerate(rows): i idx[r[0]] j idx[r[1]] for c in range(3): B[3*k c, 3*i c] -1.0 # 起点负 B[3*k c, 3*j c] 1.0 # 终点正 return B符号约定必须和常数项 l 的计算方式严格一致。如果 l 定义成观测减近似那么起点就是 −1、终点 1如果反过来定义成近似减观测符号要全部对调否则解出来的坐标改正数方向是反的最终坐标会错得离谱。另一个常见错误是索引错位B 矩阵的第 3k、3k1、3k2 行分别对应 X、Y、Z如果列方向上也用同样的排列3i、3i1、3i2一一对应即可但一旦有人在某处把点号和分量的顺序写混结果就是 X 分量的观测被算到了 Y 的未知数上。3.3 秩亏与基准约束的两种处理方式自由网的 N 矩阵秩亏 3直接求逆会报错或者给出垃圾解。两种处理方式我都用过第一种是用最小二乘的通用求解器拿最小范数解dx, *_ np.linalg.lstsq(B, l, rcondNone)它的数学含义是求满足观测方程且改正数平方和最小的那个解同时范数最小等价于固定网的重心。做内部检核很方便缺点是基准的物理含义不直观。第二种是显式施加约束构造扩展法方程# C 是约束矩阵例如固定 0 号点的三个坐标 C np.zeros((3, 3 * n_pts)) C[0, 0] C[1, 1] C[2, 2] 1.0 N B.T P B W B.T P l A np.block([[N, C.T], [C, np.zeros((3, 3))]]) sol np.linalg.solve(A, np.concatenate([W, np.zeros(3)])) dx sol[:3 * n_pts]这种做法的好处是约束的物理意义明确我确实要固定这个点而且精度评定可以直接从 N 的逆里取。坏处是如果 C 的构造有问题法方程可能依然奇异numpy 会抛 LinAlgError这时候要回头检查 C 的秩是不是 3。3.4 精度评定单位权中误差、点位误差、卡方检验解完之后不做精度评定等于白做。三个必看的量V B dx - l sigma0 np.sqrt(V.T P V / (3 * len(rows) - (3 * n_pts - 3))) Qxx np.linalg.inv(N) point_std [sigma0 * np.sqrt(Qxx[3*i, 3*i] Qxx[3*i1, 3*i1] Qxx[3*i2, 3*i2]) for i in range(n_pts)]单位权中误差 σ₀是整体拟合质量的指标。它接近 1说明协方差阵的量级和实际符合程度匹配远大于 1说明观测精度没有协方差阵声称的那么好或者附近有粗差远小于 1说明协方差阵偏保守或者人为加了比例因子。它大于 1 是常态不必惊慌。点位中误差由 Qxx 的对角块给出。注意 Qxx 是加了约束之后的逆自由网里讨论点位精度本身意义不大因为那取决于你怎么定基准。卡方检验用来判断模型是否成立VᵀPV / σ₀² 应服从自由度等于多余观测数的 χ² 分布。落在两侧的临界值之外就要考虑随机模型是否设置错了。三个指标里我最看重的是改正数的分布。把所有基线的 V 按分量画出来看如果只是随机散布在零附近、没有超过 3σ 的心里就稳了如果出现某几条边的改正数同号且偏大通常是粗差或者这些边和其他边不是同一批数据解出来的比如混了不同星历、不同时段的高度角设置。4. 从 WGS-84 落到地方坐标约束平差与转换参数的顺序问题4.1 三维无约束到二维约束中间不能跳步完整的规范流程是三步走我见过太多人直接从三维基线跳到平面坐标结果一堆莫名其妙的问题第一步三维无约束平差。目的是检核网本身的内部一致性、剔除粗差、评估相对精度。这一步的成果是相对于网重心的坐标不能直接拿去用。第二步三维约束平差。选部分已知点参与约束把网的基准对齐到目标坐标系比如 CGCS2000 或地方坐标系的空间直角坐标。这一步的关键产出是已知点兼容性检验——比较已知点坐标和约束后坐标的残差看哪些点不兼容。不兼容的点要么剔除要么降权。第三步二维约束平差。在平面坐标系里做得到最终的地方平面坐标。这一步会引入投影和方位必须先把已知点的高斯平面坐标反算到椭球面或者直接在三维空间里做平面约束。跳过第二步的风险在于你把所有已知点一股脑约束进去其中某个点有 5 cm 的误差整个网会被它拖着转而你在平面坐标里看不出任何异常。4.2 七参数与四参数选错就差一个数量级布尔沙七参数模型三个平移、三个旋转、一个尺度共七个参数。适合大范围、跨投影带、三维坐标之间的转换比如 WGS-84 到 CGCS2000 或到地方独立坐标系。平面四参数两个平移、一个旋转、一个尺度共四个参数。适合小范围、同一投影带内、二维平面坐标之间的转换比如两套城市坐标之间的对接。选择依据就两条控制点的空间分布范围和坐标系基准是否一致。控制点范围超过一个投影带大约 3° 经差或者跨省用七参数四参数在大范围里会把投影变形吃进旋转参数里越远偏得越狠。控制点集中在一个城市内、都在同一带四参数足够而且四参数需要的最少公共点更少。高差大、点位分布跨越较大高程带必须用七参数因为四参数在平面上没法处理高程方向的影响。参数的个数直接决定了最少公共点数七参数至少 3 个四参数至少 2 个。但这只是数学上的最低要求实际使用中公共点数量应该是参数个数的两倍以上否则参数解没有冗余、无法检验一个点有粗差就能把整套参数带偏。解出参数后一定要留出几个点做检核看残差是否在厘米级。4.3 边长反算与投影改化这两项不做网就永远对不上这是最容易被跳过、后果也最严重的一步。GPS 基线是从空间得到的斜距或椭球面距离要拿去和地面网比较或者输入平面坐标必须做两项归算第一项高程归算归算到参考椭球面。近似公式是 ΔS₁ ≈ −Δh²/(2S)其中 Δh 是相对于参考椭球面的高差。举个数平均高程 100 m边长 1 km改正量约 −5 mm/km。听着不大但如果测区平均高程 500 m改正量就变成 −125 mm/km也就是 1 km 差 12.5 cm这已经远超多数工程网的允许误差了。第二项高斯投影改化归算到高斯平面。近似公式是 ΔS₂ ≈ S·y²/(2R²)y 是测区到中央子午线的平均横坐标R 是地球平均曲率半径。在 y 50 km 处投影变形约 1/32000也就是 1 km 差 3 cmy 100 km 处1 km 差 13 cm。两项的符号通常相反会有一部分抵消但绝不能指望它们自己抵完。判断测区是否要做投影改化的经验值是相对变形的绝对值超过 1/40000就必须处理。碰到投影变形过大的测区常见的工程解法是建立地方独立坐标系——选一个合适的抵偿面或者中央子午线把变形控制在可接受范围内。一个真实的返工案例测区平均高程 800 m、距中央子午线 60 km两项改正叠加后边长相对误差约 1/25000。当时没做归算直接拿 GPS 坐标算出来的边和地面网差了几十厘米第一反应以为是 GPS 数据有问题查了两天闭合环、换了三次星历最后发现是根本没做投影改化。5. 大地高、正常高以及把点画到地图上的那些偏移5.1 大地高不是海拔差的是高程异常GPS 直接解算出来的是大地高 H也就是点到参考椭球面的距离。我们平时说的海拔是正常高 h是点到似大地水准面的距离。两者之间的关系是H h ζζ 就是高程异常。它的量级在全国范围内从几米到几十米不等跟地形和地壳结构有关。如果你直接把 GPS 解出的大地高当成海拔填进成果表误差就是几十米这个量级——这个错误不致命因为太好发现了但它会让后续所有基于高程的处理都错位。5.2 GPS 高程拟合能用到什么范围把 GPS 和水准结合用少量重合点拟合出测区的高程异常曲面从而把大量 GPS 点的大地高换算成正常高这是常规做法。有用的经验是二次曲面拟合至少需要 6 个均匀分布的重合点。低于这个数曲面系数解不出来或者极不稳定。拟合范围要覆盖所有待求点不要外推。外推超过重合点包络范围 10% 以外精度下降得非常快。地形起伏大的地方山地区曲面拟合基本没法用。这时改用高程异常模型比如 EGM2008 之类的重力场模型给出一个初值精度大概在分米级对一般工程够用。平原地区、重合点分布均匀、点距 2~5 km二次曲面拟合的精度通常能到 2~3 cm这是比较乐观的情况。还有一个细节水准点的高程本身可能有沉降或者历史遗留的误差。拟合前先做一次重合点之间的高程异常变化趋势检查看有没有明显的异常点别把脏数据喂进拟合。5.3 BLH 到高德坐标GCJ-02 不是一个平移量数据要落到地图上就会碰到坐标系的问题。目前常见的有三套坐标系使用场景与 WGS-84 的关系WGS-84GPS 原始输出、国际标准基准GCJ-02高德、腾讯等国内地图非线性偏移偏移量随位置变化BD-09百度地图在 GCJ-02 基础上再做一次偏移关键点是 GCJ-02 的偏移不是简单加一个常数而是一组基于经纬度的非线性变换。所以我用某个点算出来的偏移量直接加在所有点上这种做法在几公里范围内大概能用跨几十公里就完全不准了。公开实现的思路大致是这样import math a 6378245.0 ee 0.00669342162296594323 def _tf_lat(lng, lat): ret -100.0 2.0*lng 3.0*lat 0.2*lat*lat 0.1*lng*lat 0.2*math.sqrt(abs(lng)) ret (20.0*math.sin(6.0*lng*math.pi) 20.0*math.sin(2.0*lng*math.pi)) * 2.0/3.0 ret (20.0*math.sin(lat*math.pi) 40.0*math.sin(lat/3.0*math.pi)) * 2.0/3.0 ret (160.0*math.sin(lat/12.0*math.pi) 320.0*math.sin(lat*math.pi/30.0)) * 2.0/3.0 return ret def _tf_lng(lng, lat): ret 300.0 lng 2.0*lat 0.1*lng*lng 0.1*lng*lat 0.1*math.sqrt(abs(lng)) ret (20.0*math.sin(6.0*lng*math.pi) 20.0*math.sin(2.0*lng*math.pi)) * 2.0/3.0 ret (20.0*math.sin(lng*math.pi) 40.0*math.sin(lng/3.0*math.pi)) * 2.0/3.0 ret (150.0*math.sin(lng/12.0*math.pi) 300.0*math.sin(lng/30.0*math.pi)) * 2.0/3.0 return ret def wgs84_to_gcj02(lng, lat): dlat _tf_lat(lng - 105.0, lat - 35.0) dlng _tf_lng(lng - 105.0, lat - 35.0) radlat lat / 180.0 * math.pi magic math.sin(radlat) magic 1 - ee * magic * magic sqrtmagic math.sqrt(magic) dlat (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi) dlng (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi) return lng dlng, lat dlat反方向GCJ-02 回 WGS-84没有解析解通常用迭代逼近或者一次近似。这部分有两个坑值得说一是高德的接口参数顺序是经度,纬度跟平时说经纬度的顺序是反的写错一次你会看到点跑到了另一个省二是这个转换本身有 1~2 m 量级的误差需要亚米级精度的成果绝对不能靠它得用带明确参数的坐标转换。成果导出到地图上做可视化时我习惯先转成 GeoJSON 或者带表头的 CSV把点号、转换后的经纬度、点位中误差一起输出。点位中误差直接作为属性字段带上在地图上按大小渲染颜色一眼就能看出网的精度分布——比看表格快得多。6. 复盘几个真踩过的坑6.1 单位和比例因子mm²、m²、微弧度单位问题我在前面提过一次但值得再强调因为它出现的频率实在太高。三种具体表现协方差阵用 mm² 而坐标用 m角度相关的量用微弧度而代码里按弧度算比例 ppm 和纯数字混用。这三种错误的共同特征是结果看起来合理但精度指标荒谬——坐标本身对不上是运气对上了是巧合。我现在的做法是写一个单位检查函数在读入数据后立刻验证某条已知边的中误差算出来应该在 5 mm 到 2 cm 之间如果算出来是 0.000005 或者 5000就说明单位有问题。这个检查花不了十分钟能省掉一整天。6.2 已知点兼容性检验不做等于把误差焊死约束平差里最危险的操作是把所有已知点都约束上去。已知点数量多的时候看起来约束越强、结果越可靠实际上是把每个已知点的误差都硬塞进了网形里。经典的判断指标是无约束平差坐标和已知坐标之间的差值如果某个已知点的差值明显大于其他点、且方向与周围点不一致它大概率有问题。处理的顺序是先用全部已知点做一次约束平差看残差分布把残差最大的点剔除重做重复直到剩下的点残差都在限差内。剔除的点不是数据坏了可能只是它属于另一套基准、或者经历了沉降但如果强行用最终成果的责任在你。一个具体的限差参考已知点残差超过 3 倍的点位中误差或者超过规范规定的已知点限差就应该被怀疑。不同等级控制网的具体数值差别很大别照搬。6.3 观测数据本身的可信度别指望平差兜底这是我最想强调的一点也是很多人忽略的平差是数据处理的最后一道关口它只能处理符合模型假设的误差。如果观测数据本身出了问题——信号受到干扰、载噪比在某个时段异常下降、伪距残差呈现整体偏移、多路径在特定方位角上持续存在——平差出来的坐标依然会是一个数学上完美的解单位权中误差可能还特别漂亮但它和真实位置的关系已经断了。所以在平差之前我会做几件事看每个站每个时段的载噪比时间序列异常下降或者分布异常的都拎出来看伪距和相位残差的空间分布如果残差在某个方位角上系统偏大多半是多路径那个方向的观测权重应该降低把同一时段的不同解算策略跑一遍做交叉验证结果差异大就说明数据本身不稳定。另外平差结果还可以用来反查数据问题如果改正数的空间分布呈现明显的系统性方向偏移比如所有基线在某个方向上改正数都是正的那多半不是随机误差而是数据处理链路上游的某个环节出了系统性问题。这种时候回头看原始观测比在平差环节调参数有用得多。说到底GPS 网平差教会我的不是那几行最小二乘公式而是先怀疑数据再怀疑模型最后才怀疑算法这个顺序。我见过太多人在算法里绕圈子其实问题一开始就在天线高、点名或者单位上。
