搞GIS和大地测量的朋友，十有八九被“高程基准”折磨过。同一座山，RTK测出来的高程和SRTM栅格读出来的高程差十几米，LiDAR点云和LiDAR生成的DEM再比一次，又差一截。多数人第一反应是设备精度不行，其实很多时候是数据“参考面”根本就不是一回事。EGM96模型就是这个参考面中绕不开的一个名字，大量老牌DEM产品都明确声明自己以EGM96大地水准面为垂直基准。这篇帖子我打算把“用EGM96模型校正DEM数据”这件事，从原理、判断规则，到本地命令行和在线工具，完整拆开讲一遍。内容比较适合刚接触高程基准转换的新手，也适合手里攒了一堆高程数据、不知道该不该统一基准的老测绘。
1. 为什么要校正：先搞清楚三种“高”之间的差
1.1 椭球高、正高和大地水准面差距
核心公式就一个：h = H + N。
这里h是椭球高，指的是从WGS84椭球面沿法线量到地面的距离。GPS/RTK直接测出来的一般就是椭球高，这也是很多GNSS设备“不设置任何改正参数”时给出的默认高度。H是正高，也叫海拔高，是从大地水准面沿重力方向量到地面的距离。水准测量、传统测绘里的“海拔”基本就是正高这个概念。N则是大地水准面差距，也叫大地水准面高、高程异常（严格来说高程异常和正高有细微区别，这里先不细抠）。
你可以把WGS84椭球想象成一个数学规则、表面光滑的标准圆球，而大地水准面是“考虑了地球质量分布不均之后，静止海平面延伸穿过大陆的一个不规则曲面”。它不像椭球面那样能用一条公式描述，只能通过重力场模型来逼近。N就是这两个面之间在某个点上的垂直距离，不同地方差很多，不是一个小常数。
数量级上，全球范围N大概在-106米到+85米之间。中国大陆大部分地区N为负值，东部沿海可能只有负几米，到了青藏高原西缘能到负四五十米。也就是说，如果你的坐标在西部山区，直接把椭球高当成海拔用，误差可能就是几十米。这个量级对淹没分析、坡度坡向计算、土方量估算、InSAR形变监测来说，都是致命的。
1.2 EGM96模型是什么，为什么现在还在用
EGM96全称Earth Gravity Model 1996，是上世纪九十年代发布的全球地球重力场模型。它用球谐函数展开到360阶次，空间分辨率大概在0.5度（约55公里）量级，工程里常用的是15弧分网格（大约27公里采样间隔）的版本，也就是大家常见的egm96_15.gtx这类文件。
虽然这些年EGM2008（2159阶）、EGM2020（2190阶）都出来了，EGM96在行业里依然无处不在。原因很简单：大量存量数据是以它为垂直基准的。SRTM、NASADEM、ASTER GDEM、ALOS AW3D30，这些全球公开DEM产品，文档里十有八九都写着“reference to EGM96 geoid”。你如果手里拿着这些数据去做校正，反而会把本来正确的值改错。EGM96地位更像是一个“历史基准”，新老数据对它都有依赖，做数据融合时你绕不开它。
1.3 哪些DEM需要校正，哪些不需要
判断核心就一条：看产品元数据里“垂直基准”写的是什么。
我按常见数据源整理了一个速查表：
| 数据源 | 常见垂直基准 | 是否需要做EGM96校正 |
|---|---|---|
| SRTM（1弧秒/3弧秒） | EGM96 | 不需要，直接使用 |
| NASADEM | EGM96 | 不需要 |
| ASTER GDEM v3 | EGM96 | 不需要 |
| ALOS AW3D30 | EGM96 | 不需要 |
| Copernicus DEM（GLO-30/GLO-90） | EGM2008 | 要转EGM96时需转换 |
| 原始LiDAR点云（RTK/PPK解算） | 常见为WGS84椭球高 | 需要校正 |
| InSAR处理后的高程产品 | 常见为椭球高 | 需要校正 |
我见过不少朋友拿到Copernicus DEM直接拿来做测区背景，没注意它是EGM2008基准，结果和其他EGM96数据一并使用后，局部区域出现一两米的系统性偏差。所以工作流第一步，永远是确认你的数据基准，而不是急着跑算法。
1.4 工具选型：在线查询、网格下载与本地处理
校正过程中你一定会用到“查N值”和“批量改栅格”两个操作。查N值适合用在线工具，简单快；批量校正适合用GDAL或Python脚本，稳定可复现。我常用的工具如下：
| 工具 | 类型 | 用途 | 入口/查找方式 |
|---|---|---|---|
| ICGEM Calculation Service | 在线 | 全球任意点EGM96大地水准面差距计算、网格下载 | icgem.gfz-potsdam.de |
| NOAA Geoid Height Calculator | 在线 | 单点大地水准面高度速查 | geodesy.noaa.gov 搜“Geoid Height Calculator” |
| EarthScope Geoid Height Calculator | 在线 | 支持EGM96的单点查询，界面简洁 | unavco.org 搜“Geoid Height Calculator” |
| PROJ/proj-data自带EGM96网格 | 本地 | GDAL和pyproj做基准转换时自动调用 | 装proj-data包或从机构官网下载 |
| GeographicLib GeoidEval | 本地命令行 | 批量查N值、离线计算 | geographiclib.sourceforge.io |
在线工具的东西通常会改版，域名和按钮位置都会微调，所以我在下文里给的更多是“怎么找、怎么用”的方法，而不是死记一个截图。真正要长期复用的校正流程，还是要落到本地脚本上。
2. 实操前置：数据体检、网格准备与环境搭建
2.1 拿到DEM先“体检”，不要急着校
我拿到任何DEM的第一步，是一行命令看元数据：
gdalinfo dem_原始数据.tif重点看三个地方。第一，Coordinate System是什么；第二，像元大小单位是度还是米；第三，有没有NoData Value以及数值是多少。顺便用gdalinfo -stats看一眼最大最小值，心里有个数。
举个真实例子：有一次我拿到一份LiDAR生成的DEM，从文件名看像是某测区的“海拔高程”，但打开元数据发现坐标系是EPSG:4979（WGS84三维地理坐标，高程轴默认就是椭球高），而且坐标本身是用RTK直接采的。这种情况不做校正，直接拿去和在线的SRTM比较，必然对不上。很多项目报告就是“由于高程基准不一致，本次对比误差较大”草草结束，说白了就是没做这一步。
2.2 准备EGM96大地水准面网格
校正DEM的本质，是给每一个像元都配一个N值，然后把这个N从椭球高里减掉。所以你必须先有一张覆盖工作区的EGM96大地水准面网格。
最省事的办法，是让GDAL/PROJ环境里自带egm96_15.gtx。安装OSGeo4W、QGIS或通过conda装gdal时，proj-data通常会把egm96_15.gtx一起带进来。Linux下常见路径是：
/usr/share/proj/egm96_15.gtxWindows下OSGeo4W里一般是：
C:\OSGeo4W\share\proj\egm96_15.gtx如果你的环境里没有，就从ICGEM或NOAA官网下载EGM96的网格文件。下载后放到一个固定目录，例如D:\geoid\egm96_15.gtx，后续命令里写绝对路径最稳妥。
另外，ICGEM还支持直接下载GeoTIFF或文本格式的EGM96网格，你可以在网页里选择区域和分辨率，生成覆盖你测区的局部网格。这个对只做小范围项目的人来说很好用，生成的文件可以直接丢进QGIS里看N值变化。
2.3 搭一个干净的本地处理环境
本地处理我推荐两条路线。如果只是偶尔校正一次，用QGIS自带的“栅格计算器”加GDAL命令行就够了；如果经常处理，建议用Python环境。
Python环境装起来也不复杂：
pip install rasterio numpy pyproj这里特别提醒：pyproj依赖底层的PROJ库，网格文件必须能被PROJ找到。如果你发现调用垂直转换时报“Cannot find geoid grid”之类的错，多半是PROJ_LIB环境变量没有指向proj-data所在目录。可以先用下面这段代码测试网格是否可用：
from pyproj import Transformer transformer = Transformer.from_crs("EPSG:4979", "EPSG:5773", always_xy=True) lon, lat, h = 100.0, 35.0, 3000.0 lon, lat, H = transformer.transform(lon, lat, h) print(f"椭球高 h = {h}") print(f"校正后 H = {H:.3f}")代码里的EPSG:4979是WGS84三维地理坐标（椭球高），EPSG:5773是WGS84 + EGM96高程（正高）。如果你的环境正常，打印出来的H应该比h小几十米，具体取决于该点N值。
3. 完整实操：把椭球高的DEM校正到EGM96正高
3.1 方法一：重采样网格加栅格减法，直观且最好排查
这是我最推荐给普通用户的方法：把EGM96网格重采样到与DEM完全一致的范围、分辨率、投影，然后用栅格计算器做减法。
假设输入文件叫dem_ellip.tif，坐标系是WGS84地理坐标（EPSG:4326），像元大小是0.00027777778度（约30米），范围在96E到107E、26N到34N之间。先把EGM96网格裁剪重采样：
gdalwarp -overwrite \ -s_srs "EPSG:4326" -t_srs "EPSG:4326" \ -te 96 26 107 34 \ -tr 0.00027777778 0.00027777778 \ -r bilinear \ /path/to/egm96_15.gtx \ geoid_egm96.tif这条命令里，-te后面的四个数字是从DEM里读出的范围，-tr是像元大小。-r bilinear用双线性插值，因为DEM是连续面，双线性比最邻近法平滑。
做完后再执行栅格减法：
gdal_calc.py \ -A dem_ellip.tif \ -B geoid_egm96.tif \ --outfile=dem_egm96.tif \ --calc="A-B" \ --NoDataValue=-9999 \ --type=Float32 \ --overwrite为什么要A-B而不是A+B？回到公式h = H + N，现在已知椭球高h和N，求正高H，所以H = h - N。后面我会专门讲这个最容易搞反的地方。--type=Float32也建议加上，否则遇到整型DEM，减完之后的小数会被直接抹掉，表面看着没问题，实际精度已经坏了。
校正完成后立刻验证。用gdalinfo -stats看新文件的最小值、最大值和平均值，再对比原文件。如果测区N为负，那么新值应该普遍变大；如果N为正，新值应该变小。数值变化方向和量级，要和你用在线工具查到的N基本一致。
3.2 方法二：GDAL直接用geoidgrids做垂直基准变换
如果你用的是GDAL 3.x和较新版本PROJ，还有一种更“顺滑”的做法：在gdalwarp里同时指定输入和输出的坐标系，让PROJ自己完成水平投影和垂直基准的变换。
gdalwarp -overwrite \ -s_srs "EPSG:4979" \ -t_srs "EPSG:4326+5773" \ -r bilinear \ dem_ellip.tif \ dem_egm96_via_gdalwarp.tif思路是：输入是带椭球高的WGS84三维坐标（EPSG:4979），输出是WGS84二维加EGM96正高的复合坐标系（EPSG:4326+5773）。GDAL在转换过程中会检测到“高程参考面”的变化，自动去PROJ的网格目录里找EGM96网格，然后完成加减N的运算。
这个方法的好处是命令短、不容易把范围搞错；缺点是一旦环境里的网格文件缺失或者PROJ版本太老，就直接报错。如果你只是想快速试一下，可以跑跑看；但如果是生产环境，我更推荐方法一，因为每一步都看得见、查得着，出了问题能定位到具体环节。
3.3 方法三：Python脚本批量处理
当你有几十个分幅DEM要批量校正时，命令行一条条跑会崩溃。这时候可以用Python配合rasterio写一个批处理脚本。核心思路和方法一完全一样，只是把“重采样”和“减法”两步合并进代码里。
import numpy as np import rasterio from rasterio.warp import reproject, Resampling from pyproj import Transformer input_path = r"D:\data\dem_ellip.tif" geoid_path = r"D:\geoid\egm96_15.tif" output_path = r"D:\data\dem_egm96.tif" with rasterio.open(input_path) as src: dem = src.read(1) profile = src.profile.copy() dst_crs = src.crs resolution = src.res bounds = src.bounds with rasterio.open(geoid_path) as geoid_src: geoid = np.zeros((src.height, src.width), dtype=np.float32) reproject( source=geoid_src.read(1), destination=geoid, src_transform=geoid_src.transform, src_crs=geoid_src.crs, dst_transform=src.transform, dst_crs=dst_crs, resampling=Resampling.bilinear, ) dem = dem.astype(np.float32) dem[dem == profile["nodata"]] = np.nan corrected = dem - geoid profile.update(dtype="float32", nodata=-9999, compress="deflate") with rasterio.open(output_path, "w", **profile) as dst: dst.write(corrected, 1)脚本里用np.nan代替NoData参加计算，避免出现“NoData值减N”这种垃圾结果。最后输出时再统一替换成-9999并启用DEFLATE压缩，文件体积能小不少。批量处理时用os.listdir或glob遍历目录就行，我自己习惯把每个文件对应的命令信息打出来，方便写日志回溯。
3.4 用在线工具做单点校验
无论是方法一还是方法三，最后都要用已知的点验证一下结果。我一般会在测区里选三到五个点，代表山区、平原、河流附近等不同地形，用在线工具算出这些点的N值作为“标准答案”。
以ICGEM为例：打开ICGEM的Calculation Service页面，选择EGM96模型，输出量选“Geoid undulation”，输入经纬度，提交后就能得到对应的N。这个方法尤其适合验证山区，因为山区N值变化快，如果你重采样的范围或分辨率弄错了，在线工具和栅格里的值一对比立刻就露馅。
对比规则是：在线工具算的N，应该和你在geoid_egm96.tif里采样的N基本一致。两者差在10厘米以内通常算正常，如果差到几十厘米以上，建议检查你用的网格版本、重采样参数，以及在线工具是否默认用了更高阶展开。
4. 常见问题与排查技巧实录
4.1 最容易搞反：加N还是减N
这个问题出现过太多次，每次都有同事问。再帮大家理一遍。
椭球高、正高和N的关系是：
h = H + N所以：
椭球高 → 正高：H = h - N 正高 → 椭球高：h = H + N用中国大陆的例子来感受：假如某点N为-35米，RTK测出椭球高是5213米。你要得到正高，直接拿5213减掉(-35)，实际得到5248米。如果方向搞反，变成5178米，直接差70米，整个项目就废了。
我建议每个项目开始前，先找一个已知高程点的控制点，用在线工具查N，然后手工算一遍“椭球高转正高”，确认方向没问题后再跑全图。这个动作只花两分钟，能拦住绝大多数低级错误。
4.2 GDAL报“Cannot find grid”之类的错
典型报错信息里一般有egm96_15.gtx或proj.db的关键词。原因就是PROJ找不到EGM96网格。解决办法有三个：
一是确认网格文件确实存在，环境变量PROJ_LIB有没有指向proj的data目录。二是如果你用的是命令行，可以把网格路径写进命令里，比如：
-s_srs "+proj=longlat +datum=WGS84 +geoidgrids=/path/to/egm96_15.gtx"三是如果你在Python里用pyproj，装完pyproj后可以执行pyproj.datadir.get_data_dir()看它到底在找哪个目录，然后把网格文件丢进去。
这个坑比较烦，很多人以为是自己命令写错了，其实只是环境问题。
4.3 Copernicus DEM是EGM2008，怎么转成EGM96
Copernicus DEM的垂直基准在公开文档里写得很清楚，是EGM2008。如果你非要把EGM2008的高程数据转到EGM96上，可以用一个间接公式：
H96 = H2008 + N2008 - N96也就是说，需要分别算出该点EGM2008的大地水准面差距和EGM96的大地水准面差距，然后对高程做一个微调。ICGEM支持同时输出不同模型的N值，很方便。操作时先在ICGEM下载EGM2008和EGM96两张网格，分别重采样到DEM的范围，然后用栅格计算器：
--calc="A + B - C"其中A是原DEM，B是EGM2008网格，C是EGM96网格。
注意，EGM96和EGM2008在全球大部分地区差异在1米以内，但部分高山区域可能到3到5米。如果你的项目精度要求不高，可能可以忽略；但如果涉及水库库容、建筑标高这类精细计算，最好转一下。
4.4 重采样后边缘出现大片NoData
这个问题几乎每次裁剪网格都会遇到。原因很简单：原始EGM96网格的范围是全球的，而你的DEM范围是局部的，gdalwarp在裁剪时由于边界对齐问题，可能导致边缘有一圈像元落到了原始网格有效范围外。
解决办法有两个。第一个是在gdalwarp里加-wo SAMPLE_STEPS=1或者加大-te范围，让裁剪框稍微比DEM范围大一圈。第二个更简单，做完减法后，用原DEM的NoData掩膜把边缘清一遍：
gdal_calc.py -A dem_egm96.tif -B dem_原始.tif --calc="where(B==-9999, -9999, A)"在Python里也可以用np.where做类似操作。总之原则是：不要让校正过程给原始有效区域添加新的NoData。
4.5 在线工具查出来的N和本地栅格对不上
这不一定是你错了。ICGEM、NOAA这些在线工具使用的模型展开阶次、插值算法、网格版本可能和本地egm96_15.gtx有细微差别。这些差别在绝大多数区域小于10厘米，不会影响DEM校正。但如果你的测区是高山峡谷，N值空间变化又猛，那么“在线工具单点结果”和“15分网格插值结果”之间出现几十厘米偏差也正常。
如果严格要求一致，最好的办法就是保持一个数据源：下载网格时顺手把版本、分辨率、日期记录下来，在线验证时也选择相同模型版本，然后允许一定的容差。别再遇到“怎么和网上不一样”就怀疑自己算错了。
5. 收尾：一点经验
校正DEM这件事，技术上不难，难的是判断“到底要不要校、往哪个方向校”。我之前处理一个无人机LiDAR项目，项目方直接给了RTK测出来的椭球高做淹没分析，测区是一片平原，N大约是-12米，结果整个模拟水面线比真事低了12米，下游该淹的没淹，不该淹的全泡了。后来做完校正，结果完全变了个样。也是从那之后我才养成习惯：不管是下载的DEM还是团队内部生产的点云，开工第一件事永远是gdalinfo、看文档、查垂直基准，绝不信文件名里的“海拔”两个字。
最后分享一个小技巧：把你的EGM96网格文件固定存放在一个目录里，比如D:\geoid\或/data/geoid/，然后把校正DEM的几行命令稍微封装成一个小脚本。以后再来新数据，改个输入路径就能跑。这套东西用不了多少时间，却能让你在数据融合、成果发布、内外业对接时少操很多心。建议现在就打开一份DEM试试，哪怕只是查一个点的N值，也比收藏夹里吃灰强。
