企业官网建设流程全解析

简介：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网格文件，本质上是从同一组系数导出的成品，适合拿来做比对验证。文件每一行记录的字段如下表所示。

所谓「完全归一化」，是指球谐基函数 P̄nm(cosθ) 的归一化条件为 ∫₀^π P̄nm² sinθ dθ = 2(2-δ₀m)，即 m=0 时积分为 4，m>0 时积分为 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', encoding='utf-8', errors='ignore') 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]保存正弦系数。(nmax+1)*(nmax+2)//2是三角系数个数公式，对应 360 阶完整展开应该是 65281 个；如果打印结果与这个数字不一致，说明源文件被截断或混入了非标准行。类型检查放在try/except里，是为了兜住部分发行版在行尾附加额外描述文字的情况，这类尾巴用len(parts)判断拦不住，只有数值解析异常才兜得住。

加载完成后，还建议顺手做一次 C20 异常系数替换。WGS84 椭球的完全归一化 C20 大约是 -4.84169454e-4，EGM96 给的是完整引力位系数；如果要算扰动位和扰动重力异常，需要把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 接近零却并非严格为零，直接用n>=2起步求和即可。

重力异常 Δ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)(2n+1))/((n-m)(n+m))) b̄ = sqrt(((2n+1)(n+m-1)(n-m-1))/((2n-3)(n+m)(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*n+1.0) / ((n-m)*(n+m))) b_bar = np.sqrt((2.0*n+1.0)*(n+m-1.0)*(n-m-1.0) / ((2.0*n-3.0)*(n+m)*(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 倍，对百万点规模的任务，差距会从几十分钟拉到几小时。

阶数不是越高越好。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̄ₙₘ - (n+m) 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 - (n+m)P[n-1,m])/sinθ递推，n=m 时直接取极限形式。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 再跑一遍，两条结果之差就是高阶贡献量级，这也是判断是否该引入区域改正模型的依据。

本文还有配套的精品资源，点击获取