简介EGM-1996-all.rar 是一套基于 EGM96 地球重力场模型的实用计算工具面向测绘、大地测量与地球物理领域的科研人员、工程师及高校相关专业学生。利用该资源可快速求解指定位置的重力异常、高程异常、垂线偏差等关键参数无需从零手写复杂球谐展开算法直接运行程序或调用工程源码即可得到结果大幅提升数据处理与教学演示效率。压缩包共 36 个文件总体积约 1.78MB内含 C# 完整工程源码解决方案、窗体与核心类、可直接运行的 exe 可执行文件、egm96.gfc 球谐系数模型文件、演示效果图gif以及海洋重力异常等说明文档docx/txt覆盖从模型加载、坐标输入、参数计算到成果导出的完整流程。代码中已包含常用坐标和预设算例方便读者对照验证也可按需修改观测点坐标、模型阶数等参数移植到自己的项目中。目前已有 678 人学习下载是学习重力场模型应用与二次开发的实用参考资料。1. 从 EGM-1996-all.rar 说起EGM96 重力异常与垂线偏差的同一套算法EGM96 是 1996 年发布的全地球重力场模型球谐展开到 360 阶次空间分辨率约 55 公里。虽然 EGM2008、XGM2019 这些新模型已经把阶数推到 2000 阶以上但 EGM96 因为系数公开、文件格式固定、位系数与 WGS84 椭球参数天然配套仍然是重力异常改算、垂线偏差估计和高程异常内插这三类任务里被反复使用的公共参考。大多数工程人员拿到EGM-1996-all.rar之后真正需要的不是解压而是把里面的位系数读对、算准。下面就直接用这整套位系数把重力异常和垂线偏差两条计算路径完整走一遍附带可以改参数就能跑的代码。2. EGM-1996-all.rar 系数文件怎么读格式、归一化与 Python 解析2.1 EGM96 位系数格式、表头与完全归一化解压 EGM-1996-all.rar 后你会看到若干份文件。最常见的核心文件是EGM96_to360.ascii这类位系数表里面是完整的球谐展开系数其余像geoid、gravity网格文件本质上是从同一组系数导出的成品适合拿来做比对验证。文件每一行记录的字段如下表所示。字段含义单位/说明第 1 列n球谐阶数 degree无量纲整数从 0 开始第 2 列m球谐次数 order0 ≤ m ≤ n第 3 列Cnm余弦项系数完全归一化无量纲第 4 列Snm正弦项系数完全归一化无量纲第 5-6 列对应 Cnm、Snm 的标准差部分文件缺省解析时可选所谓「完全归一化」是指球谐基函数 P̄nm(cosθ) 的归一化条件为 ∫₀^π P̄nm² sinθ dθ 2(2-δ₀m)即 m0 时积分为 4m0 时积分为 2。凡是勒让德递推计算只要递推系数与这一归一化一致Cnm、Snm 不需要额外换算要防止的是拿旧版未归一化系数套归一化递推那会得到整体偏小或偏大的结果俗称「半归一化」混用错误。EGM96 官方文件标的就是完全归一化。另一个常被忽略的动作是跳过表头。有的发行版在文件开头带#注释行有的用一行说明记录成果版本解析时按字符判断过滤即可。做工程时我建议在前端加一个小peek函数打印前 10 行的字段长度分布如果某一行字段数少于 4直接跳过防止整个数组因错位解出负阶数。2.2 Python 快速读取系数到 numpy 数组读取逻辑分三步逐行过滤注释与空行、解析前四列、按C[n,m]形式填充二维数组。下面是完整实现。import numpy as np def load_egm96(fname): 读取 EGM96 位系数文件返回 C、S 二维数组及最大阶数 nmax data [] with open(fname, r, encodingutf-8, errorsignore) as f: for line in f: line line.strip() # 跳过注释行、空行与字段数不足的行 if not line or line.startswith(#): continue parts line.split() if len(parts) 4: continue try: n, m int(parts[0]), int(parts[1]) cnm, snm float(parts[2]), float(parts[3]) except ValueError: continue if m 0 or m n: continue data.append((n, m, cnm, snm)) nmax max(d[0] for d in data) C np.zeros((nmax 1, nmax 1)) S np.zeros((nmax 1, nmax 1)) for n, m, cnm, snm in data: C[n, m] cnm S[n, m] snm return C, S, nmax C, S, nmax load_egm96(EGM96_to360.ascii) print(nmax , nmax, 系数个数 , (nmax 1) * (nmax 2) // 2)该函数用numpy初始化两个(361, 361)的二维数组C[n,m]保存余弦系数S[n,m]保存正弦系数。(nmax1)*(nmax2)//2是三角系数个数公式对应 360 阶完整展开应该是 65281 个如果打印结果与这个数字不一致说明源文件被截断或混入了非标准行。类型检查放在try/except里是为了兜住部分发行版在行尾附加额外描述文字的情况这类尾巴用len(parts)判断拦不住只有数值解析异常才兜得住。加载完成后还建议顺手做一次 C20 异常系数替换。WGS84 椭球的完全归一化 C20 大约是 -4.84169454e-4EGM96 给的是完整引力位系数如果要算扰动位和扰动重力异常需要把C[2,0]替换成两者之差。这一步在下一章的重力异常计算里体现垂线偏差的球近似计算可以暂不替换因为正常场的水平分量在球近似下为零。3. 用 EGM96 球谐综合计算重力异常公式、代码与截断阶数3.1 重力异常的球谐展开式与正常场扣除计算重力异常要从扰动位 T 出发它是真实地球引力位 V 与正常椭球引力位 U 之差。在球坐标 (r, θ, λ) 下EGM96 位系数给出的扰动位展开式为T(r,θ,λ) (GM/r) · Σₙ (a/r)ⁿ · Σₘ [ΔC̄ₙₘ cos(mλ) S̄ₙₘ sin(mλ)] · P̄ₙₘ(cosθ)这里的 ΔC̄ₙₘ 表示已扣除正常椭球位系数后的异常系数θ 是地心余纬λ 是地心经度a 是参考椭球长半轴。一阶项在展开中通常被跳过因为当坐标原点取在地心时一阶项对应的质量中心偏差严格为零实际读入的 EGM96 文件里 C10、C11、S11 接近零却并非严格为零直接用n2起步求和即可。重力异常 Δg 是扰动位沿径向求导后与 2T/r 的组合代入展开式后得到Δg (GM/r²) · Σₙ (n-1)(a/r)ⁿ · Σₘ [ΔC̄ₙₘ cos(mλ) S̄ₙₘ sin(mλ)] · P̄ₙₘ(cosθ)系数 (n-1) 是理解阶数贡献的关键n 越大(a/r)ⁿ 越小所以 360 阶展开对近地面点的主要贡献来自中低阶到了卫星轨道高度高阶项迅速衰减这也是卫星重力反演只能恢复中低阶的物理原因。工程上如果测区范围小于 10°×10°把截断阶数从 360 降到 100 阶分辨率损失对区域趋势影响很小计算量却可以降一个量级。正常场扣除是最容易出错的地方。EGM96 官方位系数包含地球全部质量分布给出的 C20它与 WGS84 椭球的 C20 不是一回事需要先替换成异常系数。只做这一步替换是因为 WGS84 把正常重力场定义到 C20 这一阶更高阶带谐项在正常场里被定义为零所以不需要扣。如果跳过这一步得到的重力异常里会多出地球扁率级的系统差量级可达数百 mGal直接毁掉对比精度。3.2 连带勒让德递推与重力异常计算代码P̄ₙₘ(cosθ) 的计算采用行递推从 (n-1,m) 和 (n-2,m) 推出 (n,m)。对完全归一化系数递推常数写成ā sqrt(((2n-1)(2n1))/((n-m)(nm))) b̄ sqrt(((2n1)(nm-1)(n-m-1))/((2n-3)(nm)(n-m)))递推形式是 P̄ₙₘ ā cosθ P̄ₙ₋₁,ₘ - b̄ P̄ₙ₋₂,ₘ。对角项 P̄ₙₙ 和次对角项 P̄ₙ,ₙ₋₁ 需要先算出来再进入 m 从 0 到 n-2 的普通项循环。def legendre_row(nmax, theta): 完全归一化连带勒让德递推返回 P[n][m] P np.zeros((nmax 1, nmax 1)) ct, st np.cos(theta), np.sin(theta) P[0, 0] 1.0 if nmax 1: P[1, 0] np.sqrt(3.0) * ct P[1, 1] np.sqrt(3.0) * st for n in range(2, nmax 1): P[n, n] np.sqrt((2.0 * n 1) / (2.0 * n)) * st * P[n-1, n-1] P[n, n-1] np.sqrt(2.0 * n 1) * ct * P[n-1, n-1] for m in range(0, n - 1): a_bar np.sqrt((2.0*n-1.0)*(2.0*n1.0) / ((n-m)*(nm))) b_bar np.sqrt((2.0*n1.0)*(nm-1.0)*(n-m-1.0) / ((2.0*n-3.0)*(nm)*(n-m))) P[n, m] a_bar * ct * P[n-1, m] - b_bar * P[n-2, m] return P def gravity_anomaly_egm96(lat_deg, lon_deg, h_m, C, S, nmax): EGM96 扰动重力异常单位 mGal GM, a 3.986004415e14, 6378136.3 e2 0.00669437999014 phi, lam np.deg2rad(lat_deg), np.deg2rad(lon_deg) # 经纬高转地心直角坐标再转地心余纬 N a / np.sqrt(1 - e2 * np.sin(phi)**2) x (N h_m) * np.cos(phi) * np.cos(lam) y (N h_m) * np.cos(phi) * np.sin(lam) z (N * (1 - e2) h_m) * np.sin(phi) r np.sqrt(x*x y*y z*z) theta np.arccos(z / r) # 替换 WGS84 正常椭球场的 C20 C_use C.copy() C_use[2, 0] - -4.841694537e-4 P legendre_row(nmax, theta) m_arr np.arange(nmax 1) cos_mlam np.cos(m_arr * lam) sin_mlam np.sin(m_arr * lam) s 0.0 for n in range(2, nmax 1): row 0.0 for m in range(0, n 1): coeff C_use[n, m] * cos_mlam[m] S[n, m] * sin_mlam[m] row coeff * P[n, m] s (n - 1) * (a / r)**n * row dg GM / (r * r) * s * 1e5 # 正常重力椭球面未加高度改正 gamma 9.7803253359 * (1 0.00193185265241 * np.sin(phi)**2) / np.sqrt(1 - e2 * np.sin(phi)**2) return dg, gamma代码中主干是双层循环叠加球谐综合。内层循环对同一阶的 m 求和外层循环累加 (n-1)(a/r)ⁿ 的贡献递推结果直接作为 P̄ₙₘ 参与乘积。C_use[2, 0] - -4.841694537e-4这一行做的是「EGM96 的完整 C20」减去「WGS84 正常椭球的 C20」减完后 C_use 才是扰动位系数。dg是模型侧重力异常gamma是椭球面正常重力严格讲地面点的正常重力还要加自由空气改正但在几十米高程以下这个误差在亚 mGal 级。如果测区高差超过 1000 米建议在对比实测时显式加上正常重力高度改正。3.3 截断阶数怎么选分辨率与区域应用用 EGM96 做区域重力异常时截断阶数一般按目标分辨率定。360 阶对应约 55 公里半波长如果只需要构造尺度趋势120 阶已经足够用于局部重力勘探或航空重力处理时往往还需要更高阶的局部模型EGM96 本身就不太够。另一个判断依据是计算点密度单点 360 阶的计算量大约是 120 阶的 9 倍对百万点规模的任务差距会从几十分钟拉到几小时。截断阶数近似分辨率典型用途120约 160 km区域趋势、高程异常大尺度改正180约 110 km省级重力异常场、垂线偏差概算360约 55 km全球网格产品对照、1°×1° 重力归算阶数不是越高越好。EGM96 的高阶系数精度随阶数下降360 阶处的误差往往比 200 阶处大一个量级如果没有实测重力点约束用 300 阶以上反而会引入噪声。实践中我常把全阶计算和截断到 180 阶各跑一遍两者差异能直观反映高频段贡献量级也能帮助判断是否需要引入局部改正模型。提示判断截断阶数是否合适最直接的办法是把 360 阶与 180 阶结果都算一遍取差值均方根。若该值已经小于应用阈值说明高频段贡献可忽略后续直接用 180 阶即可。4. EGM96 垂线偏差计算ξ、η 分量的递推与符号陷阱4.1 垂线偏差公式与球坐标下的导数项垂线偏差定义为真实重力方向与正常重力方向之间的夹角。按天文大地测量习惯分解为子午圈分量 ξ 和卯酉圈分量 η单位常用角秒。由扰动位 T 求垂线偏差的球近似公式是ξ (1/(rγ)) · ∂T/∂φ η (1/(rγ cosφ)) · ∂T/∂λ注意第一个导数是关于地理纬度 φ而非余纬 θ。代入 φπ/2-θ 的关系后得到 ξ -(1/(rγ)) · ∂T/∂θ。实际编码中很多框架直接使用余纬递推所以负号被压进 ξ 的表达式。η 的表达式没有负号因为它对 λ 直接求偏导λ 增大即向东正东分量取正。推导并不复杂。把 T 的球谐展开对 θ 求导需要 dP̄ₙₘ/dθ 的递推对 λ 求导则是对 sin(mλ) 与 cos(mλ) 求导后再乘 m。两类导数项在代码里独立计算最后各自做球谐综合。4.2 勒让德导数递推与 ξ、η 的稳定计算P̄ₙₘ 关于余纬的导数递推式在 θ 不接近 0 或 π 时用dP̄ₙₘ/dθ (n cosθ P̄ₙₘ - (nm) P̄ₙ₋₁,ₘ) / sinθ这个式子来源直接实现简单但在极点附近 sinθ 趋于零数值误差会被放大。对绝对纬度大于 80° 的极区测点建议降低最大阶数到 120 阶或者改用极区专用梯度递推避免高频分量在分母上造成数值爆炸。下面是完整代码它在legendre_row之上又做了一次导数递推。def deflection_egm96(lat_deg, lon_deg, h_m, C, S, nmax): EGM96 垂线偏差返回 (xi, eta) 单位角秒 GM, a 3.986004415e14, 6378136.3 e2 0.00669437999014 phi, lam np.deg2rad(lat_deg), np.deg2rad(lon_deg) N a / np.sqrt(1 - e2 * np.sin(phi)**2) x (N h_m) * np.cos(phi) * np.cos(lam) y (N h_m) * np.cos(phi) * np.sin(lam) z (N * (1 - e2) h_m) * np.sin(phi) r np.sqrt(x*x y*y z*z) theta np.arccos(z / r) st, ct np.sin(theta), np.cos(theta) P legendre_row(nmax, theta) # 对余纬求导接近极点时做置零简化处理 dP np.zeros_like(P) for n in range(1, nmax 1): for m in range(0, n 1): if np.abs(st) 1e-10: dP[n, m] 0.0 elif n m: dP[n, m] n * ct / st * P[n, m] else: dP[n, m] (n * ct * P[n, m] - (n m) * P[n-1, m]) / st gamma 9.7803253359 * (1 0.00193185265241 * np.sin(phi)**2) / np.sqrt(1 - e2 * np.sin(phi)**2) m_arr np.arange(nmax 1) cos_mlam np.cos(m_arr * lam) sin_mlam np.sin(m_arr * lam) sum_xi 0.0 sum_eta 0.0 for n in range(2, nmax 1): row_xi 0.0 row_eta 0.0 for m in range(0, n 1): c_term C[n, m] * cos_mlam[m] S[n, m] * sin_mlam[m] d_term -C[n, m] * sin_mlam[m] S[n, m] * cos_mlam[m] row_xi c_term * dP[n, m] row_eta m * d_term * P[n, m] sum_xi (a / r)**n * row_xi sum_eta (a / r)**n * row_eta rad2arcsec 180.0 / np.pi * 3600.0 xi -GM / (r * r * gamma) * sum_xi * rad2arcsec eta GM / (r * r * gamma * st) * sum_eta * rad2arcsec return xi, eta这段代码和重力异常函数共享同一个legendre_row所以两者的球谐综合很难出现系统性不一致。dP 数组用(n cosθ P - (nm)P[n-1,m])/sinθ递推nm 时直接取极限形式。eta 的分母保留st因为球坐标里经度方向的弧长因子是 r sinθ如果这里漏掉st高纬度区 η 值会整体偏大且偏差随纬度增大而增大这是垂线偏差代码里很典型的错误。xi前面的负号来自余纬导数和纬度导数的关系不是随意加的。如果把 θ 换成地理纬度 φ 再算公式变成正值但多数球谐库的输入是余纬所以「负号加余纬」是工程默认配置。最稳妥的做法是拿已知测区标定中国东部平坦地区 ξ 一般在 -5″ 到 5″ 之间η 在 -10″ 到 10″ 之间如果一个点算出来整体符号相反优先检查负号。4.3 符号约定的交付规范垂线偏差在应用中带符号参与大量计算比如天文经纬度归算、惯导系统重力扰动补偿符号一旦传错会直接抵消正确信号。我的习惯是在交付数据时附一个 CSV 字段说明注明字段名、正方向、单位、计算模型、截断阶数和基准椭球。至少要把「ξ 正值指向北η 正值指向东」写清楚。EGM96 在全球大多数地区的垂线偏差模型值不超过 30″与实测天文大地垂线偏差之差通常能到 2-5″这已经是全球模型目前的实用边界。5. EGM96 结果验证三种对照方法与三个实操坑计算做完必须验收我习惯先做三种对照。第一种是官方网格比对EGM96 同源发布过 15′ 网格的重力异常和垂线偏差文件在测点做最近邻或双线性插值与球谐计算结果比较正常应该在 0.1 mGal 和 0.1″ 量级若出现数百 mGal 的系统差第一嫌疑是 C20 正常场没扣干净。第二种是和 EGM2008 交叉验证两边都截断到 360 阶全球大部分地区重力异常差在 1 mGal 内、垂线偏差差在 1″ 内差异远超这个范围优先检查归一化是否混用。第三种是实测点残差残差如果随高程线性走就是正常重力高度改正缺失残差随经纬度做长波变化多半是截断阶数过低造成信号泄漏。实操上三个坑最常踩。一是表头不过滤EGM-1996-all.rar 里的 ASCII 文件有的带版本说明行np.loadtxt直接读会失败或错位宁可先扫描再解析。二是极区导数递推发散纬度超过 80° 时公式分母接近零简单方案是限阶到 180用噪声换稳定。三是数组索引转置如果按C[m][n]而不是C[n][m]存储结果「大趋势对、细节错」验收时特意做一个低阶展开比对就能暴露。python check_egm96.py --file EGM96_to360.ascii \ --points 0,0 0,90 35,120 -30,60 \ --repr 360 \ --tol 0.05把--repr 360换成 180 再跑一遍两条结果之差就是高阶贡献量级这也是判断是否该引入区域改正模型的依据。本文还有配套的精品资源点击获取
EGM96球谐展开实战：重力异常与垂线偏差的完整计算实现
简介EGM-1996-all.rar 是一套基于 EGM96 地球重力场模型的实用计算工具面向测绘、大地测量与地球物理领域的科研人员、工程师及高校相关专业学生。利用该资源可快速求解指定位置的重力异常、高程异常、垂线偏差等关键参数无需从零手写复杂球谐展开算法直接运行程序或调用工程源码即可得到结果大幅提升数据处理与教学演示效率。压缩包共 36 个文件总体积约 1.78MB内含 C# 完整工程源码解决方案、窗体与核心类、可直接运行的 exe 可执行文件、egm96.gfc 球谐系数模型文件、演示效果图gif以及海洋重力异常等说明文档docx/txt覆盖从模型加载、坐标输入、参数计算到成果导出的完整流程。代码中已包含常用坐标和预设算例方便读者对照验证也可按需修改观测点坐标、模型阶数等参数移植到自己的项目中。目前已有 678 人学习下载是学习重力场模型应用与二次开发的实用参考资料。1. 从 EGM-1996-all.rar 说起EGM96 重力异常与垂线偏差的同一套算法EGM96 是 1996 年发布的全地球重力场模型球谐展开到 360 阶次空间分辨率约 55 公里。虽然 EGM2008、XGM2019 这些新模型已经把阶数推到 2000 阶以上但 EGM96 因为系数公开、文件格式固定、位系数与 WGS84 椭球参数天然配套仍然是重力异常改算、垂线偏差估计和高程异常内插这三类任务里被反复使用的公共参考。大多数工程人员拿到EGM-1996-all.rar之后真正需要的不是解压而是把里面的位系数读对、算准。下面就直接用这整套位系数把重力异常和垂线偏差两条计算路径完整走一遍附带可以改参数就能跑的代码。2. EGM-1996-all.rar 系数文件怎么读格式、归一化与 Python 解析2.1 EGM96 位系数格式、表头与完全归一化解压 EGM-1996-all.rar 后你会看到若干份文件。最常见的核心文件是EGM96_to360.ascii这类位系数表里面是完整的球谐展开系数其余像geoid、gravity网格文件本质上是从同一组系数导出的成品适合拿来做比对验证。文件每一行记录的字段如下表所示。字段含义单位/说明第 1 列n球谐阶数 degree无量纲整数从 0 开始第 2 列m球谐次数 order0 ≤ m ≤ n第 3 列Cnm余弦项系数完全归一化无量纲第 4 列Snm正弦项系数完全归一化无量纲第 5-6 列对应 Cnm、Snm 的标准差部分文件缺省解析时可选所谓「完全归一化」是指球谐基函数 P̄nm(cosθ) 的归一化条件为 ∫₀^π P̄nm² sinθ dθ 2(2-δ₀m)即 m0 时积分为 4m0 时积分为 2。凡是勒让德递推计算只要递推系数与这一归一化一致Cnm、Snm 不需要额外换算要防止的是拿旧版未归一化系数套归一化递推那会得到整体偏小或偏大的结果俗称「半归一化」混用错误。EGM96 官方文件标的就是完全归一化。另一个常被忽略的动作是跳过表头。有的发行版在文件开头带#注释行有的用一行说明记录成果版本解析时按字符判断过滤即可。做工程时我建议在前端加一个小peek函数打印前 10 行的字段长度分布如果某一行字段数少于 4直接跳过防止整个数组因错位解出负阶数。2.2 Python 快速读取系数到 numpy 数组读取逻辑分三步逐行过滤注释与空行、解析前四列、按C[n,m]形式填充二维数组。下面是完整实现。import numpy as np def load_egm96(fname): 读取 EGM96 位系数文件返回 C、S 二维数组及最大阶数 nmax data [] with open(fname, r, encodingutf-8, errorsignore) as f: for line in f: line line.strip() # 跳过注释行、空行与字段数不足的行 if not line or line.startswith(#): continue parts line.split() if len(parts) 4: continue try: n, m int(parts[0]), int(parts[1]) cnm, snm float(parts[2]), float(parts[3]) except ValueError: continue if m 0 or m n: continue data.append((n, m, cnm, snm)) nmax max(d[0] for d in data) C np.zeros((nmax 1, nmax 1)) S np.zeros((nmax 1, nmax 1)) for n, m, cnm, snm in data: C[n, m] cnm S[n, m] snm return C, S, nmax C, S, nmax load_egm96(EGM96_to360.ascii) print(nmax , nmax, 系数个数 , (nmax 1) * (nmax 2) // 2)该函数用numpy初始化两个(361, 361)的二维数组C[n,m]保存余弦系数S[n,m]保存正弦系数。(nmax1)*(nmax2)//2是三角系数个数公式对应 360 阶完整展开应该是 65281 个如果打印结果与这个数字不一致说明源文件被截断或混入了非标准行。类型检查放在try/except里是为了兜住部分发行版在行尾附加额外描述文字的情况这类尾巴用len(parts)判断拦不住只有数值解析异常才兜得住。加载完成后还建议顺手做一次 C20 异常系数替换。WGS84 椭球的完全归一化 C20 大约是 -4.84169454e-4EGM96 给的是完整引力位系数如果要算扰动位和扰动重力异常需要把C[2,0]替换成两者之差。这一步在下一章的重力异常计算里体现垂线偏差的球近似计算可以暂不替换因为正常场的水平分量在球近似下为零。3. 用 EGM96 球谐综合计算重力异常公式、代码与截断阶数3.1 重力异常的球谐展开式与正常场扣除计算重力异常要从扰动位 T 出发它是真实地球引力位 V 与正常椭球引力位 U 之差。在球坐标 (r, θ, λ) 下EGM96 位系数给出的扰动位展开式为T(r,θ,λ) (GM/r) · Σₙ (a/r)ⁿ · Σₘ [ΔC̄ₙₘ cos(mλ) S̄ₙₘ sin(mλ)] · P̄ₙₘ(cosθ)这里的 ΔC̄ₙₘ 表示已扣除正常椭球位系数后的异常系数θ 是地心余纬λ 是地心经度a 是参考椭球长半轴。一阶项在展开中通常被跳过因为当坐标原点取在地心时一阶项对应的质量中心偏差严格为零实际读入的 EGM96 文件里 C10、C11、S11 接近零却并非严格为零直接用n2起步求和即可。重力异常 Δg 是扰动位沿径向求导后与 2T/r 的组合代入展开式后得到Δg (GM/r²) · Σₙ (n-1)(a/r)ⁿ · Σₘ [ΔC̄ₙₘ cos(mλ) S̄ₙₘ sin(mλ)] · P̄ₙₘ(cosθ)系数 (n-1) 是理解阶数贡献的关键n 越大(a/r)ⁿ 越小所以 360 阶展开对近地面点的主要贡献来自中低阶到了卫星轨道高度高阶项迅速衰减这也是卫星重力反演只能恢复中低阶的物理原因。工程上如果测区范围小于 10°×10°把截断阶数从 360 降到 100 阶分辨率损失对区域趋势影响很小计算量却可以降一个量级。正常场扣除是最容易出错的地方。EGM96 官方位系数包含地球全部质量分布给出的 C20它与 WGS84 椭球的 C20 不是一回事需要先替换成异常系数。只做这一步替换是因为 WGS84 把正常重力场定义到 C20 这一阶更高阶带谐项在正常场里被定义为零所以不需要扣。如果跳过这一步得到的重力异常里会多出地球扁率级的系统差量级可达数百 mGal直接毁掉对比精度。3.2 连带勒让德递推与重力异常计算代码P̄ₙₘ(cosθ) 的计算采用行递推从 (n-1,m) 和 (n-2,m) 推出 (n,m)。对完全归一化系数递推常数写成ā sqrt(((2n-1)(2n1))/((n-m)(nm))) b̄ sqrt(((2n1)(nm-1)(n-m-1))/((2n-3)(nm)(n-m)))递推形式是 P̄ₙₘ ā cosθ P̄ₙ₋₁,ₘ - b̄ P̄ₙ₋₂,ₘ。对角项 P̄ₙₙ 和次对角项 P̄ₙ,ₙ₋₁ 需要先算出来再进入 m 从 0 到 n-2 的普通项循环。def legendre_row(nmax, theta): 完全归一化连带勒让德递推返回 P[n][m] P np.zeros((nmax 1, nmax 1)) ct, st np.cos(theta), np.sin(theta) P[0, 0] 1.0 if nmax 1: P[1, 0] np.sqrt(3.0) * ct P[1, 1] np.sqrt(3.0) * st for n in range(2, nmax 1): P[n, n] np.sqrt((2.0 * n 1) / (2.0 * n)) * st * P[n-1, n-1] P[n, n-1] np.sqrt(2.0 * n 1) * ct * P[n-1, n-1] for m in range(0, n - 1): a_bar np.sqrt((2.0*n-1.0)*(2.0*n1.0) / ((n-m)*(nm))) b_bar np.sqrt((2.0*n1.0)*(nm-1.0)*(n-m-1.0) / ((2.0*n-3.0)*(nm)*(n-m))) P[n, m] a_bar * ct * P[n-1, m] - b_bar * P[n-2, m] return P def gravity_anomaly_egm96(lat_deg, lon_deg, h_m, C, S, nmax): EGM96 扰动重力异常单位 mGal GM, a 3.986004415e14, 6378136.3 e2 0.00669437999014 phi, lam np.deg2rad(lat_deg), np.deg2rad(lon_deg) # 经纬高转地心直角坐标再转地心余纬 N a / np.sqrt(1 - e2 * np.sin(phi)**2) x (N h_m) * np.cos(phi) * np.cos(lam) y (N h_m) * np.cos(phi) * np.sin(lam) z (N * (1 - e2) h_m) * np.sin(phi) r np.sqrt(x*x y*y z*z) theta np.arccos(z / r) # 替换 WGS84 正常椭球场的 C20 C_use C.copy() C_use[2, 0] - -4.841694537e-4 P legendre_row(nmax, theta) m_arr np.arange(nmax 1) cos_mlam np.cos(m_arr * lam) sin_mlam np.sin(m_arr * lam) s 0.0 for n in range(2, nmax 1): row 0.0 for m in range(0, n 1): coeff C_use[n, m] * cos_mlam[m] S[n, m] * sin_mlam[m] row coeff * P[n, m] s (n - 1) * (a / r)**n * row dg GM / (r * r) * s * 1e5 # 正常重力椭球面未加高度改正 gamma 9.7803253359 * (1 0.00193185265241 * np.sin(phi)**2) / np.sqrt(1 - e2 * np.sin(phi)**2) return dg, gamma代码中主干是双层循环叠加球谐综合。内层循环对同一阶的 m 求和外层循环累加 (n-1)(a/r)ⁿ 的贡献递推结果直接作为 P̄ₙₘ 参与乘积。C_use[2, 0] - -4.841694537e-4这一行做的是「EGM96 的完整 C20」减去「WGS84 正常椭球的 C20」减完后 C_use 才是扰动位系数。dg是模型侧重力异常gamma是椭球面正常重力严格讲地面点的正常重力还要加自由空气改正但在几十米高程以下这个误差在亚 mGal 级。如果测区高差超过 1000 米建议在对比实测时显式加上正常重力高度改正。3.3 截断阶数怎么选分辨率与区域应用用 EGM96 做区域重力异常时截断阶数一般按目标分辨率定。360 阶对应约 55 公里半波长如果只需要构造尺度趋势120 阶已经足够用于局部重力勘探或航空重力处理时往往还需要更高阶的局部模型EGM96 本身就不太够。另一个判断依据是计算点密度单点 360 阶的计算量大约是 120 阶的 9 倍对百万点规模的任务差距会从几十分钟拉到几小时。截断阶数近似分辨率典型用途120约 160 km区域趋势、高程异常大尺度改正180约 110 km省级重力异常场、垂线偏差概算360约 55 km全球网格产品对照、1°×1° 重力归算阶数不是越高越好。EGM96 的高阶系数精度随阶数下降360 阶处的误差往往比 200 阶处大一个量级如果没有实测重力点约束用 300 阶以上反而会引入噪声。实践中我常把全阶计算和截断到 180 阶各跑一遍两者差异能直观反映高频段贡献量级也能帮助判断是否需要引入局部改正模型。提示判断截断阶数是否合适最直接的办法是把 360 阶与 180 阶结果都算一遍取差值均方根。若该值已经小于应用阈值说明高频段贡献可忽略后续直接用 180 阶即可。4. EGM96 垂线偏差计算ξ、η 分量的递推与符号陷阱4.1 垂线偏差公式与球坐标下的导数项垂线偏差定义为真实重力方向与正常重力方向之间的夹角。按天文大地测量习惯分解为子午圈分量 ξ 和卯酉圈分量 η单位常用角秒。由扰动位 T 求垂线偏差的球近似公式是ξ (1/(rγ)) · ∂T/∂φ η (1/(rγ cosφ)) · ∂T/∂λ注意第一个导数是关于地理纬度 φ而非余纬 θ。代入 φπ/2-θ 的关系后得到 ξ -(1/(rγ)) · ∂T/∂θ。实际编码中很多框架直接使用余纬递推所以负号被压进 ξ 的表达式。η 的表达式没有负号因为它对 λ 直接求偏导λ 增大即向东正东分量取正。推导并不复杂。把 T 的球谐展开对 θ 求导需要 dP̄ₙₘ/dθ 的递推对 λ 求导则是对 sin(mλ) 与 cos(mλ) 求导后再乘 m。两类导数项在代码里独立计算最后各自做球谐综合。4.2 勒让德导数递推与 ξ、η 的稳定计算P̄ₙₘ 关于余纬的导数递推式在 θ 不接近 0 或 π 时用dP̄ₙₘ/dθ (n cosθ P̄ₙₘ - (nm) P̄ₙ₋₁,ₘ) / sinθ这个式子来源直接实现简单但在极点附近 sinθ 趋于零数值误差会被放大。对绝对纬度大于 80° 的极区测点建议降低最大阶数到 120 阶或者改用极区专用梯度递推避免高频分量在分母上造成数值爆炸。下面是完整代码它在legendre_row之上又做了一次导数递推。def deflection_egm96(lat_deg, lon_deg, h_m, C, S, nmax): EGM96 垂线偏差返回 (xi, eta) 单位角秒 GM, a 3.986004415e14, 6378136.3 e2 0.00669437999014 phi, lam np.deg2rad(lat_deg), np.deg2rad(lon_deg) N a / np.sqrt(1 - e2 * np.sin(phi)**2) x (N h_m) * np.cos(phi) * np.cos(lam) y (N h_m) * np.cos(phi) * np.sin(lam) z (N * (1 - e2) h_m) * np.sin(phi) r np.sqrt(x*x y*y z*z) theta np.arccos(z / r) st, ct np.sin(theta), np.cos(theta) P legendre_row(nmax, theta) # 对余纬求导接近极点时做置零简化处理 dP np.zeros_like(P) for n in range(1, nmax 1): for m in range(0, n 1): if np.abs(st) 1e-10: dP[n, m] 0.0 elif n m: dP[n, m] n * ct / st * P[n, m] else: dP[n, m] (n * ct * P[n, m] - (n m) * P[n-1, m]) / st gamma 9.7803253359 * (1 0.00193185265241 * np.sin(phi)**2) / np.sqrt(1 - e2 * np.sin(phi)**2) m_arr np.arange(nmax 1) cos_mlam np.cos(m_arr * lam) sin_mlam np.sin(m_arr * lam) sum_xi 0.0 sum_eta 0.0 for n in range(2, nmax 1): row_xi 0.0 row_eta 0.0 for m in range(0, n 1): c_term C[n, m] * cos_mlam[m] S[n, m] * sin_mlam[m] d_term -C[n, m] * sin_mlam[m] S[n, m] * cos_mlam[m] row_xi c_term * dP[n, m] row_eta m * d_term * P[n, m] sum_xi (a / r)**n * row_xi sum_eta (a / r)**n * row_eta rad2arcsec 180.0 / np.pi * 3600.0 xi -GM / (r * r * gamma) * sum_xi * rad2arcsec eta GM / (r * r * gamma * st) * sum_eta * rad2arcsec return xi, eta这段代码和重力异常函数共享同一个legendre_row所以两者的球谐综合很难出现系统性不一致。dP 数组用(n cosθ P - (nm)P[n-1,m])/sinθ递推nm 时直接取极限形式。eta 的分母保留st因为球坐标里经度方向的弧长因子是 r sinθ如果这里漏掉st高纬度区 η 值会整体偏大且偏差随纬度增大而增大这是垂线偏差代码里很典型的错误。xi前面的负号来自余纬导数和纬度导数的关系不是随意加的。如果把 θ 换成地理纬度 φ 再算公式变成正值但多数球谐库的输入是余纬所以「负号加余纬」是工程默认配置。最稳妥的做法是拿已知测区标定中国东部平坦地区 ξ 一般在 -5″ 到 5″ 之间η 在 -10″ 到 10″ 之间如果一个点算出来整体符号相反优先检查负号。4.3 符号约定的交付规范垂线偏差在应用中带符号参与大量计算比如天文经纬度归算、惯导系统重力扰动补偿符号一旦传错会直接抵消正确信号。我的习惯是在交付数据时附一个 CSV 字段说明注明字段名、正方向、单位、计算模型、截断阶数和基准椭球。至少要把「ξ 正值指向北η 正值指向东」写清楚。EGM96 在全球大多数地区的垂线偏差模型值不超过 30″与实测天文大地垂线偏差之差通常能到 2-5″这已经是全球模型目前的实用边界。5. EGM96 结果验证三种对照方法与三个实操坑计算做完必须验收我习惯先做三种对照。第一种是官方网格比对EGM96 同源发布过 15′ 网格的重力异常和垂线偏差文件在测点做最近邻或双线性插值与球谐计算结果比较正常应该在 0.1 mGal 和 0.1″ 量级若出现数百 mGal 的系统差第一嫌疑是 C20 正常场没扣干净。第二种是和 EGM2008 交叉验证两边都截断到 360 阶全球大部分地区重力异常差在 1 mGal 内、垂线偏差差在 1″ 内差异远超这个范围优先检查归一化是否混用。第三种是实测点残差残差如果随高程线性走就是正常重力高度改正缺失残差随经纬度做长波变化多半是截断阶数过低造成信号泄漏。实操上三个坑最常踩。一是表头不过滤EGM-1996-all.rar 里的 ASCII 文件有的带版本说明行np.loadtxt直接读会失败或错位宁可先扫描再解析。二是极区导数递推发散纬度超过 80° 时公式分母接近零简单方案是限阶到 180用噪声换稳定。三是数组索引转置如果按C[m][n]而不是C[n][m]存储结果「大趋势对、细节错」验收时特意做一个低阶展开比对就能暴露。python check_egm96.py --file EGM96_to360.ascii \ --points 0,0 0,90 35,120 -30,60 \ --repr 360 \ --tol 0.05把--repr 360换成 180 再跑一遍两条结果之差就是高阶贡献量级这也是判断是否该引入区域改正模型的依据。本文还有配套的精品资源点击获取
