简介：基于EGM96模型计算已知位置重力异常、高程异常与垂线偏差的实用工具包，面向大地测量、地球物理专业学生以及需要处理重力场数据的工程技术人员。资源以C#源码工程为核心，共36个文件，除项目源码外还包含可执行程序、球谐系数数据文件、坐标文本样例、运行界面图片与Word说明文档，其中调试符号、界面资源等工程辅助文件也一并保留，便于二次编译与界面定制；压缩包仅1.78MB，整体结构紧凑，方便直接运行、修改参数并嵌入到自己的项目中。已有679人学习下载。借助这份资源可完整获得Visual Studio工程、egm96.gfc系数文件、海洋重力异常分析笔记及点位坐标示例；既能快速算出指定位置的重力场参数，也能对照源码和文档理解EGM96解算流程、异常量转换关系与数据组织方式，是学习垂线偏差和重力异常解算的高质量参考样例。
1. 拿到 EGM96 重力场模型包之后，最值得做的三件事
做大地测量和物探的人，对“EGM-1996-all.rar”这类压缩包通常又爱又恨：解压出来一堆 .cof、.dat、.txt 文件，名字长得像天书，没有 README 的时候连从哪下手都不知道。这个包实际是 EGM96 地球重力场模型的完整系数与辅助数据集合，覆盖全球 360 阶球谐展开，能直接算出高程异常、垂线偏差、重力异常三大类物理量。测绘里做 GNSS 高程拟合、InSAR 形变解算时去掉长波误差，物探里做区域重力场背景场扣除，用的都是同一套底层系数。这个方向值不值得投入，取决于你手里有没有需要绝对重力基准的数据——只要涉及跨区域拼接或高程转换，EGM96 这套老模型反而比某些新模型更稳，因为它的参考椭球和 WGS84 完全一致，很多历史数据也是基于它平差的。下面从原理到复现逐步拆开讲。
2. EGM96 的系数文件到底存了什么：球谐系数、参考椭球与截断阶数
2.1 从 .cof 文件看 EGM96 的物理含义
EGM96 的数据包里，核心文件是 EGM96.cof，每一行是一个球谐系数项。格式大致是：阶数 n、次数 m、规格化余弦系数 C̄nm、规格化正弦系数 S̄nm，外加一个误差项。不要小看这几列数字，整个重力场模型就是靠它们逼近的。地球重力位 V 在球坐标下展开为无穷级数，实际计算时截断到有限阶。EGM96 展开到 360 阶，对应空间分辨率大约是半波长 55 公里（也就是 360 阶对应的最短波长约 111 公里，半波长约 55 公里）。所以用它来研究比 50 公里更小的局部重力异常是没有意义的，这一点很多初学者会忽略，拿到模型就把所有信号都当真实场，最后和实测数据一对比就翻车。
系数文件里前几行是 0 阶和 1 阶项。0 阶项代表地球总质量和 GM 值，1 阶项理论上为零——因为坐标原点取在地球质心。如果你从某个渠道下载的系数文件 1 阶项不是零，说明数据在传输或转换过程中出了问题，这是第一个可以做的快速校验。
2.2 规格化与未规格化的边界条件
球谐系数分为“完全规格化”（fully normalized）和“未规格化”（unnormalized）两种表达。EGM96 官方发布的是完全规格化系数，计算勒让德函数时也必须用完全规格化的连带勒让德函数 P̄nm。很多人把系数从文件读出来后，直接用未规格化的 Pnm 公式去算，结果后面阶数越高偏差越大。判断方法很简单：完全规格化后，所有阶数的 P̄nm 平方在球面上的积分平均值为 1，而未规格化的 Pnm 数值随阶数急剧增大或振荡。写代码时建议直接用标准递推公式算 P̄nm，不要自己从定义式积分。
EGM96 的参考椭球参数是：GM = 3986004.415×10⁸ m³/s²，a = 6378136.3 m，扁率 f = 1/298.257223563。注意这里的 a 和 WGS84 椭球的长半轴 6378137.0 m 差了 0.7 米，虽然很小，但在垂线偏差计算里会影响正常重力值，进而影响重力异常的基准。如果你想和 CGCS2000 框架下的实测数据对比，最好把椭球参数统一到目标坐标系再做差值，不要直接拿 EGM96 的绝对输出当结论。
2.3 下载包里的辅助文件怎么用
除了 .cof 系数文件，包里通常还有 geoid 格网文件（如 egm96_15min.bin 或类似格式）、垂线偏差格网、重力异常格网以及一版说明文档。格网文件是预先算好的离散值，适合快速查值；系数文件适合自己编程算任意点位。两者的关系是：格网值就是用系数文件按一定网格算出来的。如果你只需要零星几个点的高程异常，直接用系数文件算比双线性插值格网更准，尤其是在地形起伏大的区域，格网的 15 分分辨率（约 28 公里）会平滑掉局部信号。反过来，如果你要生成一张区域等值线图，格网文件算得快，系数文件逐点算会慢得多。
3. 用 Python 从球谐系数计算高程异常和重力异常：最小可运行实现
3.1 读取 EGM96.cof 的代码与字段说明
先写一个最简读取函数。EGM96.cof 每行包含 n、m、C̄nm、S̄nm 和两个误差项，用空格分隔。建议用 pandas 读，但为了在没装 pandas 的服务器上也能跑，我用标准库实现：
def load_egm96_cof(filepath): coeffs = {} with open(filepath, 'r', encoding='ascii') as f: for line in f: parts = line.split() if len(parts) < 4: continue n = int(parts[0]) m = int(parts[1]) cnm = float(parts[2]) snm = float(parts[3]) # 误差列对计算无用，忽略 coeffs[(n, m)] = (cnm, snm) return coeffs这个函数返回一个字典，键是 (n, m) 元组，值是对应系数。注意 EG M96 的系数表是按 n 从小到大、m 从 0 到 n 排列的，所以直接顺序读即可。如果文件里有注释行，要留意跳过；部分重新打包过的版本第一行可能是文件头，建议加载后检查键 (0,0) 是否存在。正规数据包里 N 最大到 360，但有时也会附一个 361 阶以上的截断版本，实际使用时应确认自己需要多大阶数——点数不够时算高阶项纯属浪费算力。
3.2 完全规格化连带勒让德函数的递推
很多实现用标准递推公式，但写错归一化系数是最常见的翻车点。下面这段代码采用数值稳定的列递推，从低阶向高阶推进：
import math def plm_bar(n_max, theta): """ 计算 sin(theta) 对应的完全规格化连带勒让德函数值。 theta 为余纬（0 到 pi），返回二维数组 P[n][m]。 """ P = [[0.0] * (n_max + 1) for _ in range(n_max + 1)] sin_theta = math.sin(theta) cos_theta = math.cos(theta) # P00 P[0][0] = 1.0 if n_max == 0: return P # P10, P11 P[1][0] = math.sqrt(3.0) * cos_theta P[1][1] = math.sqrt(3.0) * sin_theta for n in range(2, n_max + 1): # m = 0 项: 利用 Pn0 递推 P[n][0] = math.sqrt((2.0 * n - 1.0) / n) * cos_theta * P[n-1][0] \ - math.sqrt((n - 1.0) / n) * P[n-2][0] # m = n 项: 对角线 P[n][n] = math.sqrt((2.0 * n + 1.0) / (2.0 * n)) * sin_theta * P[n-1][n-1] # 内部项 m = 1..n-1 for m in range(1, n): a = math.sqrt(((2.0 * n - 1.0) * (2.0 * n + 1.0)) / ((n - m) * (n + m))) b = math.sqrt(((2.0 * n + 1.0) * (n - m - 1.0) * (n + m - 1.0)) / ((n - m) * (n + m) * (2.0 * n - 3.0))) P[n][m] = a * cos_theta * P[n-1][m] - b * P[n-2][m] return P这个递推公式是工程上最常用的版本，系数已经包含了完全规格化因子。注意输入是余纬 theta（地心余纬），不是地理纬度。地心余纬和地理纬度的差别在高精度计算里不能忽略：当你要算全球任意一点时，应该先把大地纬度转换成地心纬度，再取余纬。转换公式为 tan(phi_geocentric) = (1 - e²) * tan(phi_geodetic)，e 是第一偏心率。用错纬度造成的误差在垂线偏差分量上可达角秒量级。
3.3 计算高程异常 N 的完整函数
高程异常 N 是大地水准面到参考椭球面的距离，表达式为 N = GM/(rγ) * ΣΣ (a/r)^n * (C̄nm cos mλ + S̄nm sin mλ) * P̄nm(sin φ)，其中 γ 是正常重力值。写成 Python 如下：
R_EGM96 = 6378136.3 GM_EGM96 = 3986004.415e8 WGS84_E2 = 6.69437999014e-3 def geoid_height(lat_deg, lon_deg, coeffs, n_max=360): lat = math.radians(lat_deg) lon = math.radians(lon_deg) # 地心纬度与地心距离（近似用椭球面） sin_lat_gdl = math.sin(lat) cos_lat_gdl = math.cos(lat) N_r = R_EGM96 / math.sqrt(1 - WGS84_E2 * sin_lat_gdl**2) x = N_r * cos_lat_gdl * math.cos(lon) y = N_r * cos_lat_gdl * math.sin(lon) z = N_r * (1 - WGS84_E2) * sin_lat_gdl r = math.sqrt(x*x + y*y + z*z) # 地心余纬 theta = math.acos(z / r) # 正常重力 γ（Somigliana 公式简化） sin_phi2 = sin_lat_gdl**2 gamma = 9.7803253359 * (1 + 0.00193185265241 * sin_phi2) / \ math.sqrt(1 - WGS84_E2 * sin_phi2) P = plm_bar(n_max, theta) s = 0.0 for n in range(0, n_max + 1): factor_n = (R_EGM96 / r) ** n row_sum = 0.0 for m in range(0, n + 1): c, s_m = coeffs.get((n, m), (0.0, 0.0)) row_sum += (c * math.cos(m * lon) + s_m * math.sin(m * lon)) * P[n][m] s += factor_n * row_sum N = GM_EGM96 / (r * gamma) * s return N代码里先算出地心距 r 和地心余纬，再遍历球谐级数。注意这里没有加 0 阶项的完整形式，但 (0,0) 系数本身就包含了主要项。正常重力 γ 用的 Somigliana 公式是闭式解，比用正常椭球位展开更精确，而且和 EGM96 的椭球参数匹配得更好。如果你算的区域在低纬度，r 的差异对最终结果影响很小；但在极地，如果直接用大地纬度代替地心纬度做余纬，N 的误差最大能到几分米。这一点是在做极区 GNSS 高程转换时最容易忽略的坑。
3.4 重力异常的计算和单位陷阱
重力异常定义为实测重力值减去正常重力值，在球谐模型中用扰动位 T 对 r 求导得到。简化公式为 Δg = GM/r² * Σ (n-1) (a/r)^n * Σ (C̄nm cos mλ + S̄nm sin mλ) * P̄nm(sin φ)。注意前面多了一个 (n-1) 因子，这是扰动位求导的结果。下面这个函数直接返回 mGal（1 Gal = 1 cm/s²，1 mGal = 10⁻⁵ m/s²）：
def gravity_anomaly(lat_deg, lon_deg, coeffs, n_max=360): lat = math.radians(lat_deg) lon = math.radians(lon_deg) sin_lat_gdl = math.sin(lat) N_r = R_EGM96 / math.sqrt(1 - WGS84_E2 * sin_lat_gdl**2) x = N_r * math.cos(lat) * math.cos(lon) y = N_r * math.cos(lat) * math.sin(lon) z = N_r * (1 - WGS84_E2) * sin_lat_gdl r = math.sqrt(x*x + y*y + z*z) theta = math.acos(z / r) P = plm_bar(n_max, theta) s = 0.0 for n in range(2, n_max + 1): factor_n = (R_EGM96 / r) ** n * (n - 1) row_sum = 0.0 for m in range(0, n + 1): c, s_m = coeffs.get((n, m), (0.0, 0.0)) row_sum += (c * math.cos(m * lon) + s_m * math.sin(m * lon)) * P[n][m] s += factor_n * row_sum dg = GM_EGM96 / (r * r) * s # 转成 mGal: 1 m/s^2 = 100 Gal = 100000 mGal return dg * 1e5从 n=2 开始是因为 n=1 项在质心坐标系下为 0，n=0 项对应的是质量项，求导后为常数，通常不放进异常里。如果你发现算出来的重力异常整体偏大或偏小几百 mGal，先检查是不是把 n=0 项也加进去了。另外注意文件里的 C̄20 量级是 10⁻³，而高阶项量级是 10⁻⁶ 甚至更小，如果程序用单精度浮点变量累加，高阶项会被舍入误差吃掉，必须用双精度。
4. 垂线偏差从 EGM96 提取：分量为纲，格网与点计算两种路线
4.1 垂线偏差的物理定义与东向、北向分量
垂线偏差是重力方向与椭球法线方向之间的夹角，分为子午圈分量 ξ（南北方向）和卯酉圈分量 η（东西方向）。它和高程异常 N 的关系是：ξ = -∂N/∂φ / (M + h)，η = -∂N/∂λ / ((N_phi + h) cos φ)，其中 M 是子午圈曲率半径，N_phi 是卯酉圈曲率半径。从 EGM96 计算垂线偏差有两条路：一条是直接用球谐系数求 N 对纬度和经度的偏导数，另一条是从包里的垂线偏差格网双线性插值。前面的数学关系更准，但计算量大；后者查表快，精度受格网分辨率限制。实际做区域垂线偏差测量时，通常是先拿 EGM96 算一个“模型垂线偏差”作为背景场，再用天文大地测量或 GPS/水准得到实测值做差，所以模型本身的分辨率决定了你能提取的局部信号上限。
4.2 用差分法从高程异常格网生成垂线偏差
如果包里附带了 geoid 格网，最快的方法是差分。例如格网分辨率为 0.25 度（约 27 公里），北向分量近似为 [-N(i+1,j)+N(i-1,j)] / (2 * d_phi * M)，东向分量为 [-N(i,j+1)+N(i,j-1)] / (2 * d_lon * N_phi * cos φ)。这段逻辑写成 Python 很直接：
def deflection_from_grid(N_grid, lat_grid, lon_grid): dlat = lat_grid[1] - lat_grid[0] dlon = lon_grid[1] - lon_grid[0] xi = np.zeros_like(N_grid) eta = np.zeros_like(N_grid) for i in range(1, len(lat_grid) - 1): phi = math.radians(lat_grid[i]) M = 6335439.0 / (1 - WGS84_E2 * math.sin(phi)**2)**1.5 N_phi = R_EGM96 / math.sqrt(1 - WGS84_E2 * math.sin(phi)**2) for j in range(1, len(lon_grid) - 1): dN_dphi = (N_grid[i+1, j] - N_grid[i-1, j]) / (2 * math.radians(dlat)) dN_dlon = (N_grid[i, j+1] - N_grid[i, j-1]) / (2 * math.radians(dlon)) xi[i, j] = -dN_dphi / M eta[i, j] = -dN_dlon / (N_phi * math.cos(phi)) return xi, eta这里的 M 和 N_phi 用了近似公式，如果追求更高精度，应该用 WGS84 椭球的精确曲率半径计算公式。差分法的边界上算不了，所以最后一行一列的值是无效的。从代码写法也能看出来，北向分量 xi 和东向分量 eta 的单位是弧度，输出时需要乘以 206265 转成角秒。一个 0.25 度格网差分得到的垂线偏差，在东向分量上误差通常在 1 角秒以内，但在地形剧烈地区会偏大，因为差分相当于人为平滑了短波信号。
4.3 直接球谐求导的高精度路线
如果你不需要生成格网，只计算离散点位的高精度垂线偏差，直接对 N 的球谐表达式求偏导数更干净。公式推导不展开，关键是最后要处理对单项 (C̄nm cos mλ + S̄nm sin mλ) * P̄nm(sin φ) 的求导。对纬度求导涉及 P̄nm 对 sin φ 的导数，可以借助递推关系 dP̄nm/dφ = -[ (n+1) * P̄n+1,m - (n+m+1) * P̄n,m ] / sin φ 之类的关系，但实现时容易在分母出现 sin φ 为零的极点问题。曾有工程团队在北纬 89 度以上的测区用这个公式直接算，结果在极点附近溢出，后来改用复数形式的递推才解决。我个人的建议是：如果不是做极区或全球格网，用数值差分就够了；如果必须做极区，改用余纬 θ 作为自变量，公式里 sin θ 在 θ=0 处也为零，但在极区重力场本身变化小，可以用局部坐标系转换绕开。
下面是一个用数值偏导实现垂线偏差点计算的例子，没有极区问题——因为它在经纬网格上做局部差分，而不是求解析导数：
def deflection_point(lat_deg, lon_deg, coeffs, step=0.01): lat_r = math.radians(lat_deg) lon_r = math.radians(lon_deg) h1 = step # 角度步长，单位度 # 北向分量：N 随纬度变化 N_north = geoid_height(lat_deg + h1, lon_deg, coeffs) N_south = geoid_height(lat_deg - h1, lon_deg, coeffs) M = 6335439.0 / (1 - WGS84_E2 * math.sin(lat_r)**2)**1.5 xi = -(N_north - N_south) / (2 * math.radians(h1) * M) # 东向分量：N 随经度变化 N_east = geoid_height(lat_deg, lon_deg + h1, coeffs) N_west = geoid_height(lat_deg, lon_deg - h1, coeffs) N_phi = R_EGM96 / math.sqrt(1 - WGS84_E2 * math.sin(lat_r)**2) eta = -(N_east - N_west) / (2 * math.radians(h1) * N_phi * math.cos(lat_r)) return math.degrees(xi) * 3600, math.degrees(eta) * 3600这个方法的步长 0.01 度对应地面约 1 公里，对 360 阶模型来说已经足够密。步长太大，会平滑掉高阶层面的短波信号；步长太小，N 的浮点误差会被放大，输出噪声变大。在普通区域建议先试 0.01 度，如果输出结果出现明显锯齿，再增大到 0.05 度。运行时间方面，对单点做四次 geoid_height 计算，每次遍历约 6.5 万个系数项，在普通 PC 上耗时不到一秒，满足工程查询需求。
5. 从下载到出图的避坑清单：格式、椭球、单位与截断
5.1 坑：系数文件读出来 C̄00 对不上
现象：用包里的系数算出来的高程异常在赤道附近有 50 米左右的整体偏移，而且所有点都偏移相同量。插值格网和系数计算结果完全对不上。
原因：C̄00 值读错或漏读。EGM96 的 C̄00 是规范化的地球引力位系数，正确值接近 -1.0，但有些重打包版本把单位改成了非规格化，或者文件头有一行说明没被解析。另外，有的实现把 GM 值分开处理，公式里额外再乘 1，和直接使用系数文件里的 0 阶项不一致。
解决：在加载后打印 coeffs[(0,0)]，正常情况下应该约等于 -0.484165...（具体数值可以在说明文档里查到）。如果不是，就把代码里 geoid_height 的公式补上对应的阶项处理，或者换一份未修改的官方系数文件。这是最容易排查也最容易被忽视的问题。
5.2 坑：把 EGM96 高程异常直接当 CGCS2000 高程异常用
现象：在华北某测区，GNSS 大地高减去 EGM96 高程异常得到的正常高，与水准测量结果对比系统性差 20-40 厘米。
原因：EGM96 大地水准面是相对 WGS84 椭球的，而 CGCS2000 的高程基准是似大地水准面，两者在理论上存在系统差。更实际的原因是 EGM96 模型本身在中国区域的偏差就有分米级，因为它的全球数据同化时东亚地区的地面重力数据稀疏。这一点和模型阶数无关，换 EGM2008 也一样存在，只是量级不同。
解决：用区域内的 GPS/水准点做高程异常拟合残差修正，或者使用 EGM96 与区域似大地水准面模型的差值格网做改正。不要指望单纯提高截断阶数能消除这个系统差，它属于模型长波误差。
5.3 坑：mGal 和 10⁻⁵ m/s² 混用导致量级偏差
现象：计算的重力异常在山区普遍比实测值小一个量级，但趋势正确。
原因：代码里把重力单位换算写成了 dg * 1e2（转成 Gal）而不是 * 1e5（转成 mGal）。在打印调试时看到数字是几十，视觉上觉得“差不多”，但和实测比较就对不上。
解决：统一用 mGal 作为单位。实测重力异常通常范围在 ±100 mGal 左右，EGM96 模型输出的中纬度区域异常一般在 -50 到 +50 mGal 之间。如果看到几千的绝对值，一定换算错了。建议在函数里显式注释单位，并在输出时强制打印单位。
5.4 坑：截断阶数不一致导致边缘振荡
现象：同一个测区，用 360 阶和用 180 阶算出来的高程异常差值，在山地达到 30 厘米，而在平原只有 2 厘米。有人误以为高阶项是噪声，直接降到 50 阶用，结果区域中心和平坦区域没问题，边界区域出现波浪形。
原因：球谐模型在截断处，如果系数不是自然衰减到零，会引入“吉布斯振荡”现象。EGM96 的 360 阶系数在高阶部分量级确实变小，但不是每个 m 都单调递减，强行截断会产生等效于加了一个矩形窗的效果。
解决：不要随意降低截断阶数。如果为了计算效率想降阶，至少用 180 阶，并且只用于远场计算；近区必须用 360 阶。另外，算区域平均值时，截断误差会在边界处最明显，这时应把计算范围向外扩 0.5 度再裁剪结果。
5.5 坑：格网文件插值时用了经纬度等距双线性
现象：用包里的 geoid 格网在 60 度以上高纬度地区做双线性插值，插值结果误差高达 1 米，而低纬度只有 5 厘米。
原因：等经纬度格网在高纬度地区经线收敛，实际地面网格从接近正方形变成一个东西向被压缩的扁矩形，双线性插值的权重失真。
解决：改用球面插值，或者在高纬度将经度步长按 cos φ 修正。最简单的方式是先把经纬度换算成地心经纬度，再做双线性插值，能恢复大部分精度。
6. 验证 EGM96 计算结果的两种可靠方法：实测点对比与跨模型差值检查
拿到一套 EGM96 计算程序后，先不要急着往报告里放数字，用一个最小验证集做交叉确认。第一种验证是实测点对比。找测区内 20 个以上的 GPS/水准联合点，用 GPS 大地高减去水准正常高得到实测高程异常，与 EGM96 计算结果比较。差值的中误差如果小于 10 厘米，说明程序实现正确，模型在这个区域的表现也符合预期；如果中误差大于 30 厘米，先回去检查椭球参数和巨维项，不要怀疑是模型不好——程序错的概率远大于模型的误差。第二种验证是跨模型对比。互联网上能下载到 EGM2008 的 2190 阶系数，拿同样点位分别用 EGM96 和 EGM2008 计算高程异常，两者差值在平原应该在 20-30 厘米以内，山区可能到 1-2 米。如果差值大到 5 米以上，大概率是其中一个模型的参考框架没对齐或者读取时跳行。EGM96 的 360 阶截断信号在 EGM2008 面前是低分辨率模型，但两者在长波部分应该高度一致。把差值画成等值线图，能看到典型的短波斑点，说明两套模型的短波部分有差异，但整体趋势一致，这是正常现象。
我自己的经验是：对垂线偏差计算结果的验证最麻烦，因为很少有人手头有实测的天文大地垂线偏差数据。一个变通办法是用 EGM96 格网生成一组点位，把这些点位当“真值”，再用自己的独立系数计算程序算一遍，看差值是否接近零。如果差值在 0.05 角秒以内，说明程序和格网是一致的，程序基本可信。不要直接信任下载包里的 exe 或二进制程序，很多是多年前在 32 位系统下编译的，在新机器上可能跑出未定义行为。最后补充一个实用习惯：每次计算都输出一行元信息，包括椭球参数来源、截断阶数、系数文件 MD5 的前 8 位、计算时刻。等三个月后回看数据，你会感谢这个习惯——很多莫名其妙的结果差异，最后都是因为用了不同版本的系数文件或改了参数没记录。希望这些步骤能让你少走几趟弯路，拿到 EGM 压缩包后第一晚就能算出可信的物理量。
本文还有配套的精品资源，点击获取
