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

前阵子做高程传递项目，甲方给的水准高程和RTK测出来的椭球高差了将近三十厘米。我一开始以为是杆子没立直，重新对中整平、换基站重测，结果还是对不上。最后查了当地的大地水准面模型，才发现问题是坐标系转换时少做了“大地水准面改正”。从那天起，我把大地水准面球谐展开公式从头到尾啃了一遍，也把为什么GNSS高程不能直接当地面高程这件事彻底弄明白了。这篇文章想把公式、实现步骤和常见坑一起讲清楚，适合正在做GNSS高程转换、重力场建模或者刚接触物理大地测量学的朋友，看完可以直接上手算一个点的大地水准面高N，也搞清楚项目里那几厘米、几十厘米到底差在哪。

这里先给一个结论：大地水准面球谐展开公式，本质上就是用一组正弦、余弦和勒让德函数的加权和，去逼近真实地球重力场引起的等位面起伏。每个权重就是重力场模型里的位系数，模型阶数越高，能表达的空间尺度越小。理解了这件事，后续所有计算都顺理成章。

1. 同一个点，三个高程：为什么要绕到球谐展开

1.1 椭球高、正高、正常高，先统一说法

很多做工程测量的人手里有GNSS接收机，测出来的是WGS84椭球高。椭球高是相对于参考椭球面的高度，参考椭球就是一个规则旋转椭球，比如WGS84、GRS80。水准仪测出来的则是正高，也就是相对于大地水准面的高度。大地水准面不是规则椭球，它是不规则但处处与重力方向垂直的等位面。

所以GNSS椭球高h减去正高H，得到的就是大地水准面高N：

h = 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_{n=2}^{N_{max}}\left(\frac{a}{r}\right)^n \sum_{m=0}^{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就是点质量的引力位。真实地球因为自转和质量不均，形状不是球体，所以需要叠加高阶项。

二阶项n=2里最重要的一项是m=0的C̄_20，它对应地球的扁率，量级在10⁻⁴到10⁻³，是所有非球形项里最大的。n=3、n=4这些项就开始描述南北不对称、梨形效应等更细致的质量分布。阶数越高，描述的空间细节越小。比如n=10的项描述的是全球范围四五千公里尺度的质量异常，n=360的项则描述约一百公里尺度的信号，EGM2008这类模型最高到2190阶，对应约9公里半波长分辨率。

这就像一个地球的“傅里叶展开”，低阶项是模糊的整体轮廓，高阶项是越来越清晰的局部纹理。

2.2 带谐项、田谐项、扇谐项的几何直觉

在展开式里，P̄_nm(cosθ)乘以cos(mλ)或sin(mλ)的组合，把整个球面上的质量异常分成了不同空间图案：

• m=0：带谐项，只随纬度变化，不随经度变化。它描述的是东西方向均匀的纬向条带，比如C̄_20、C̄_40。
• 0<m<n：田谐项，在经纬度方向都有变化，是棋盘状的补丁结构。
• m=n：扇谐项，只在经度方向快速变化，像西瓜皮一样纵向分割。

每一项系数的大小决定了该图案对重力位的贡献。大地水准面高就是所有这些图案按权重叠加后的结果。所以看一份重力场模型的系数文件，本质上就是看地球质量分布在各个空间尺度上的能量分布。

2.3 阶次与空间分辨率：n=360意味着什么

球谐函数在球面上的空间波长大约为：

$$\lambda_n \approx \frac{2\pi a}{n}$$

半波长就是能分辨的最小尺度：

$$\Delta x \approx \frac{\pi a}{n}$$

比如：

所以，如果你只用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.00000000000000e+00 S 2 1 -1.47055875595157e-10 0.00000000000000e+00 C 3 0 9.57254135202790e-07 0.00000000000000e+00

这里C和S后面跟着阶数n、次数m，然后是两个浮点数，第一个是系数值，第二个是误差或标准差，实际计算只取第一个。

下载时还要确认两个东西：参考椭球和潮汐系统。EGM2008默认锚定WGS84椭球，ICGEM上新模型通常会标出推荐椭球。潮汐约定也必须看，不同模型可能用零潮、无潮，直接混用会有厘米级差异。

坐标转换方面需要注意：GNSS测出来的是大地纬度φ，而球谐展开用的是地心余纬θ。先在给定椭球下把(φ, λ, h)转成地心直角坐标(X, Y, Z)，再用：

$$r = \sqrt{X^2+Y^2+Z^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 是地心余纬，单位弧度 返回形状为 (nmax+1, nmax+1) 的二维数组 """ 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[m+1][m] if m <= nmax - 1: P[m+1, 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

注意m=1的对角线项要手动给，不能从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})(2n+1)\frac{(n-m)!}{(n+m)!}}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在线算一遍做交叉验证，再开始批量处理。另外记得把坐标转换、勒让德递推、位系数求和拆成独立函数分别测试。球谐展开公式本身不复杂，真正难的，是在大地测量、地球物理和工程坐标之间来回切换的时候，确保每一步都没搞错归一化、参考椭球和潮汐约定。这几点盯住了，几厘米的精度完全做得到。