大地水准面球谐展开公式解析：从原理到GNSS高程转换实践

前阵子做高程传递项目甲方给的水准高程和RTK测出来的椭球高差了将近三十厘米。我一开始以为是杆子没立直重新对中整平、换基站重测结果还是对不上。最后查了当地的大地水准面模型才发现问题是坐标系转换时少做了“大地水准面改正”。从那天起我把大地水准面球谐展开公式从头到尾啃了一遍也把为什么GNSS高程不能直接当地面高程这件事彻底弄明白了。这篇文章想把公式、实现步骤和常见坑一起讲清楚适合正在做GNSS高程转换、重力场建模或者刚接触物理大地测量学的朋友看完可以直接上手算一个点的大地水准面高N也搞清楚项目里那几厘米、几十厘米到底差在哪。这里先给一个结论大地水准面球谐展开公式本质上就是用一组正弦、余弦和勒让德函数的加权和去逼近真实地球重力场引起的等位面起伏。每个权重就是重力场模型里的位系数模型阶数越高能表达的空间尺度越小。理解了这件事后续所有计算都顺理成章。1. 同一个点三个高程为什么要绕到球谐展开1.1 椭球高、正高、正常高先统一说法很多做工程测量的人手里有GNSS接收机测出来的是WGS84椭球高。椭球高是相对于参考椭球面的高度参考椭球就是一个规则旋转椭球比如WGS84、GRS80。水准仪测出来的则是正高也就是相对于大地水准面的高度。大地水准面不是规则椭球它是不规则但处处与重力方向垂直的等位面。所以GNSS椭球高h减去正高H得到的就是大地水准面高Nh H N这个N就是大地水准面相对于参考椭球面的起伏全球范围大概在-110m到90m之间。如果直接把椭球高当成海拔在N比较大的区域就会像开头说的那样差出去几十厘米甚至更多。有些国家用正常高系统引入的是一个类似但不等同的“似大地水准面”对应的高程异常用ζ表示。工程上用N还是ζ取决于高程框架但背后的物理量都是靠同一个重力场模型展开来算的。1.2 大地水准面等于一个“等位面”所以要从位函数说起大地水准面的定义是与静止平均海水面重合并延伸入陆地内部的等位面。既然是等位面它的形状就取决于地球重力位W在空间里怎么分布。W分成两个部分引力位V取决于地球内部质量分布离心位地球自转引起的惯性力对应的位。在外部空间引力位V满足拉普拉斯方程。对拉普拉斯方程做分离变量得到的通解就是球谐级数。所以地球重力场位函数可以用球谐展开来表达而大地水准面高N可以通过扰动位T和正常重力γ联系起来。扰动位T定义为真实重力位W与正常重力位U的差T W - U布隆斯公式给出了T和N的关系N T / γ这就是整个“大地水准面球谐展开公式”的核心链条先获得扰动位的球谐展开再除以正常重力值得到大地水准面高。1.3 公式长什么样在球坐标下扰动位T的球谐展开形式为$$T(r,\theta,\lambda) \frac{GM}{r}\sum_{n2}^{N_{max}}\left(\frac{a}{r}\right)^n \sum_{m0}^{n}\left(\Delta\bar{C}{nm}\cos m\lambda \Delta\bar{S}{nm}\sin m\lambda\right)\bar{P}_{nm}(\cos\theta)$$其中r是计算点到地心的距离θ是地心余纬北极取0λ是经度GM是地球引力常数a是参考椭球长半轴n是阶数m是次数ΔC̄_nm、ΔS̄_nm是模型的完全归一化位系数通常已经扣除了参考椭球的正常位贡献P̄_nm(cosθ)是完全归一化缔合勒让德函数。这个公式看起来长但每一项都有物理含义。下一节我们把它拆开看。2. 公式里的每一项都在讲什么位系数就是“地球形状的数字指纹”2.1 从零阶到高阶质量、扁率与细节如果不加任何非球形项GM/r就是点质量的引力位。真实地球因为自转和质量不均形状不是球体所以需要叠加高阶项。二阶项n2里最重要的一项是m0的C̄_20它对应地球的扁率量级在10⁻⁴到10⁻³是所有非球形项里最大的。n3、n4这些项就开始描述南北不对称、梨形效应等更细致的质量分布。阶数越高描述的空间细节越小。比如n10的项描述的是全球范围四五千公里尺度的质量异常n360的项则描述约一百公里尺度的信号EGM2008这类模型最高到2190阶对应约9公里半波长分辨率。这就像一个地球的“傅里叶展开”低阶项是模糊的整体轮廓高阶项是越来越清晰的局部纹理。2.2 带谐项、田谐项、扇谐项的几何直觉在展开式里P̄_nm(cosθ)乘以cos(mλ)或sin(mλ)的组合把整个球面上的质量异常分成了不同空间图案m0带谐项只随纬度变化不随经度变化。它描述的是东西方向均匀的纬向条带比如C̄_20、C̄_40。0mn田谐项在经纬度方向都有变化是棋盘状的补丁结构。mn扇谐项只在经度方向快速变化像西瓜皮一样纵向分割。每一项系数的大小决定了该图案对重力位的贡献。大地水准面高就是所有这些图案按权重叠加后的结果。所以看一份重力场模型的系数文件本质上就是看地球质量分布在各个空间尺度上的能量分布。2.3 阶次与空间分辨率n360意味着什么球谐函数在球面上的空间波长大约为$$\lambda_n \approx \frac{2\pi a}{n}$$半波长就是能分辨的最小尺度$$\Delta x \approx \frac{\pi a}{n}$$比如阶数n全波长约(km)半波长约(km)220030100155801040053611105551802221113601115672056282160199所以如果你只用EGM2008截断到360阶算N你就丢了所有小于约56公里的重力场细节。在某些地形起伏大、质量异常集中的区域丢掉这些高频信号会造成厘米级甚至分米级的差异。2.4 绝对位系数和差分位系数最容易搞混的一步球谐展开本身是描述完整重力位V的。但计算扰动位T时要从V里减去正常椭球对应的正常位U。所以最终使用的系数有两种绝对位系数完整重力位的球谐系数差分位系数真实位系数减去参考椭球正常位系数也就是ΔC̄_nm、ΔS̄_nm。很多在线模型比如ICGEM下载服务会明确提供“difference coefficients”。如果你是直接下载完整位系数记得自己扣除参考椭球项否则算出来的T会包含一个巨大的椭球背景场N会差到米级甚至更大。使用前先看说明文件这一句话能帮你省半天排查时间。3. 从公式到数字手写一个计算N的小工具3.1 数据准备下载模型并确认坐标系统计算前先准备一份重力场模型系数。推荐用EGM2008或EIGEN-6C4去ICGEM官网下载.gfc格式文件。打开文件后每一行大概是C 2 0 -4.84169391702101e-04 0.00000000000000e00 S 2 1 -1.47055875595157e-10 0.00000000000000e00 C 3 0 9.57254135202790e-07 0.00000000000000e00这里C和S后面跟着阶数n、次数m然后是两个浮点数第一个是系数值第二个是误差或标准差实际计算只取第一个。下载时还要确认两个东西参考椭球和潮汐系统。EGM2008默认锚定WGS84椭球ICGEM上新模型通常会标出推荐椭球。潮汐约定也必须看不同模型可能用零潮、无潮直接混用会有厘米级差异。坐标转换方面需要注意GNSS测出来的是大地纬度φ而球谐展开用的是地心余纬θ。先在给定椭球下把(φ, λ, h)转成地心直角坐标(X, Y, Z)再用$$r \sqrt{X^2Y^2Z^2}$$$$\theta \arccos(Z/r)$$这样比直接用大地纬度代入更稳妥也顺便回避了地心纬度和大地纬度的转换公式。3.2 缔合勒让德函数的递推实现完全归一化缔合勒让德函数P̄_nm(x)是公式里计算量最大的部分。直接按定义算阶乘会溢出工程上都用递推。我的做法是先算到最大阶数Nmax存成一个二维数组方便后面循环。核心Python代码可以这样写import numpy as np def legendre_norm(nmax, theta): 计算完全归一化缔合勒让德函数 Pbar[n][m] theta 是地心余纬单位弧度 返回形状为 (nmax1, nmax1) 的二维数组 cos_t np.cos(theta) sin_t np.sin(theta) P np.zeros((nmax 1, nmax 1)) P[0, 0] 1.0 if nmax 1: P[1, 0] np.sqrt(3.0) * cos_t P[1, 1] np.sqrt(3.0) * sin_t for m in range(0, nmax 1): # 对角线递推 P[m][m] if m 2: P[m, m] np.sqrt((2.0 * m 1.0) / (2.0 * m)) * sin_t * P[m-1, m-1] # 次对角线递推 P[m1][m] if m nmax - 1: P[m1, m] np.sqrt(2.0 * m 3.0) * cos_t * P[m, m] # 一般递推 P[n][m] for n in range(m 2, nmax 1): a np.sqrt((4.0 * n * n - 1.0) / (n * n - m * m)) b np.sqrt(((n - 1.0) * (n - 1.0) - m * m) / (4.0 * (n - 1.0) * (n - 1.0) - 1.0)) P[n, m] a * (cos_t * P[n-1, m] - b * P[n-2, m]) return P注意m1的对角线项要手动给不能从P[0][0]用递推推出来否则会发现结果和理论值对不上。这是很多人第一次写勒让德递推时踩的坑。3.3 主循环从位系数到扰动位再到N有了P̄_nm表求和就简单了。假设你已经把模型系数读成两个二维数组C和S并且是差分系数主循环如下def geoid_height(lat_deg, lon_deg, h_ell, C, S, GM, a, nmax): # 1. 大地坐标转地心直角坐标这里用WGS84椭球为例 f 1.0 / 298.257223563 e2 f * (2.0 - f) phi np.radians(lat_deg) lam np.radians(lon_deg) N_phi a / np.sqrt(1.0 - e2 * np.sin(phi)**2) X (N_phi h_ell) * np.cos(phi) * np.cos(lam) Y (N_phi h_ell) * np.cos(phi) * np.sin(lam) Z (N_phi * (1.0 - e2) h_ell) * np.sin(phi) r np.sqrt(X*X Y*Y Z*Z) theta np.arccos(Z / r) # 2. 计算勒让德函数 P legendre_norm(nmax, theta) # 3. 计算扰动位 T T 0.0 for n in range(2, nmax 1): s 0.0 for m in range(0, n 1): angle m * lam s (C[n, m] * np.cos(angle) S[n, m] * np.sin(angle)) * P[n, m] T (a / r) ** n * s T GM / r * T # 4. 正常重力近似公式WGS84 sin_phi np.sin(phi) gamma0 9.7803253359 * (1.0 0.00193185265241 * sin_phi**2) \ / np.sqrt(1.0 - 0.00669437999014 * sin_phi**2) N T / gamma0 return N几点说明这里计算的是大地水准面高N用的是布隆斯公式T/γ如果模型给的是完整位系数而不是差分系数循环里要先把参考椭球的系数减掉h_ell是椭球高如果你只有海拔不知道椭球高可以先取0算一版N对几十米级别的高程不敏感多数情况下误差在毫米量级。3.4 用ICGEM或GMT验证你的结果代码写完不能直接信先用官方工具验算几个点。ICGEM的在线计算服务可以选模型、选坐标直接输出gravity field functionals包括geoid undulation N。挑两三个点比如经度120°、纬度30°附近对比你的结果和ICGEM结果。如果一致到毫米级代码基本没问题如果差很多优先检查是否用了差分系数纬度是不是转成了地心余纬勒让德递推有没有溢出。我自己第一次写的时候就是漏了扣参考椭球项结果N比ICGEM大了将近100米。那种错误从数值上非常明显但搜索半天不容易发现。4. 为什么别人算的和你差一截五个常见坑4.1 纬度类型用错球谐展开函数的自变量是余纬θ但这个θ是地心余纬不是测绘里常用的大地纬度。如果直接把大地纬度φ换成90°-φ代入结果在高纬度地区可能偏出几公里对N的影响能达到几十厘米到米级。原因是大地纬度和地心纬度差异最大在45°左右可到0.19°对应地面距离约20公里对中长波段的重力场求和来说相位误差非常大。正确做法还是先转地心直角坐标再用acos(Z/r)求θ。这样一步到位也顺便把r拿到手。4.2 归一化约定没对齐位系数有未归一化、部分归一化、完全归一化三种常见约定。EGM2008和ICGEM下载的模型基本都是完全归一化但很多教材和老代码用的是未归一化系数。完全归一化P̄_nm和普通缔合勒让德P_nm的关系大约是$$\bar{P}{nm} \sqrt{(2-\delta{0m})(2n1)\frac{(n-m)!}{(nm)!}}P_{nm}$$如果你用普通P_nm去乘以完全归一化系数结果会差很多反过来也一样。所以拿到别人的代码先看他有没有在勒让德函数里做归一化别盲目复用。4.3 参考椭球和潮汐系统不一致大地水准面高N是相对于某个参考椭球的。同一个重力场模型锚定WGS84和GRS80算出来的N会有差别。虽然两个椭球参数差异很小但在高精度应用里厘米级差异不可忽略。潮汐问题更容易被忽视。重力场模型可能在“零潮系统”“无潮系统”或“平均潮系统”下构建理论上它们之间的差异可以达到分米级。如果你的工程高原采用正常高系统最好选用与该系统协调一致的模型版本并在文档里写明采用的是哪种约定。4.4 球近似到底够不够公式里用了r和θ本质上是球坐标展开。但地球是个椭球严格做法应该用椭球谐函数。大多数重力场模型提供的是球谐系数所以实际计算普遍采用球近似即在公式里直接代入地心r、θ。球近似的误差有多大对中低阶项影响很小对高阶项有一定影响。实践下来计算大地水准面高的误差通常不超过几厘米而且主要在起伏大的区域。如果你要的是毫米级结果就需要引入椭球改正项或者使用专门的椭球谐展开算法。4.5 截断和地形效应所有模型都有最大阶数。截断到Nmax意味着忽略了更短波长的重力场信号。在地形陡峭、质量异常集中的区域比如高山峡谷截断误差可能很大。而且球谐展开在原理上只适用于外部无质量空间在地表以下不收敛。所以在高山区直接用球谐模型算N效果往往不如在低海拔平原地区。工程上处理这个问题的常见办法是先用球谐模型算一个长波背景场再用地形质量模型做高频改正再用局部GNSS水准数据做残差拟合。这一步到位就是区域大地水准面精化。5. 这个公式在生产里的真正用法5.1 GNSS高程转换最直接的应用就是GNSS高程转换。已知椭球高h用重力场模型算出N正常高H就是H h - N在EGM2008模型覆盖较好的地区直接用全球模型算N精度大约在十厘米到几十厘米量级取决于区域重力场复杂度。如果让模型N内插到测点再做一次简单的残差改正可以把精度推到厘米级。实际操作时不要去每个点都跑一遍球谐展开而是先用球谐展开算一个规则网格的大地水准面高比如1分或者30秒间隔再在测点做双线性内插。这样效率高精度损失也不大。5.2 区域大地水准面精化如果要进一步提高精度就需要GNSS水准控制点。每个控制点上的实测N是N_实测 h_GNSS - H_水准全球模型给出N_模型两者之差ΔN包含了模型误差、长波误差和局部高频信号。用已知控制点上的ΔN做曲面拟合比如多项式拟合或多面函数然后内插到未知点上对N_模型做改正。这个流程就是经典的“移去-恢复法”移去控制点上计算ΔN拟合用函数拟合ΔN空间趋势恢复未知点的N N_模型 拟合得到的ΔN。这样处理后区域大地水准面的精度可以做到亚厘米级。5.3 卫星重力和时变信号的延伸球谐展开不止能算静态N。GRACE、GRACE-FO和GOCE等卫星重力任务反演出的球谐系数是随时间变化的。对这些时变系数做同样的求和得到的就是大地水准面高的时间变化。通过这个方法可以监测大范围地下水储量变化、冰盖质量变化甚至洋流带来的质量迁移。我后来给客户做GNSS高程转换项目时已经不满足于只用几个控制点拟合而是把全球模型系数、区域重力数据和GPS水准一起做了联合精化。球谐展开公式始终是底层的主心骨其他所有改正项都是围着它转的。最后分享一个自己的习惯不管用哪个模型、哪份代码先找两三个已知点用ICGEM在线算一遍做交叉验证再开始批量处理。另外记得把坐标转换、勒让德递推、位系数求和拆成独立函数分别测试。球谐展开公式本身不复杂真正难的是在大地测量、地球物理和工程坐标之间来回切换的时候确保每一步都没搞错归一化、参考椭球和潮汐约定。这几点盯住了几厘米的精度完全做得到。