搞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这类文件。虽然这些年EGM20082159阶、EGM20202190阶都出来了EGM96在行业里依然无处不在。原因很简单大量存量数据是以它为垂直基准的。SRTM、NASADEM、ASTER GDEM、ALOS AW3D30这些全球公开DEM产品文档里十有八九都写着“reference to EGM96 geoid”。你如果手里拿着这些数据去做校正反而会把本来正确的值改错。EGM96地位更像是一个“历史基准”新老数据对它都有依赖做数据融合时你绕不开它。1.3 哪些DEM需要校正哪些不需要判断核心就一条看产品元数据里“垂直基准”写的是什么。我按常见数据源整理了一个速查表数据源常见垂直基准是否需要做EGM96校正SRTM1弧秒/3弧秒EGM96不需要直接使用NASADEMEGM96不需要ASTER GDEM v3EGM96不需要ALOS AW3D30EGM96不需要Copernicus DEMGLO-30/GLO-90EGM2008要转EGM96时需转换原始LiDAR点云RTK/PPK解算常见为WGS84椭球高需要校正InSAR处理后的高程产品常见为椭球高需要校正我见过不少朋友拿到Copernicus DEM直接拿来做测区背景没注意它是EGM2008基准结果和其他EGM96数据一并使用后局部区域出现一两米的系统性偏差。所以工作流第一步永远是确认你的数据基准而不是急着跑算法。1.4 工具选型在线查询、网格下载与本地处理校正过程中你一定会用到“查N值”和“批量改栅格”两个操作。查N值适合用在线工具简单快批量校正适合用GDAL或Python脚本稳定可复现。我常用的工具如下工具类型用途入口/查找方式ICGEM Calculation Service在线全球任意点EGM96大地水准面差距计算、网格下载icgem.gfz-potsdam.deNOAA 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:4979WGS84三维地理坐标高程轴默认就是椭球高而且坐标本身是用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_xyTrue) 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 \ --outfiledem_egm96.tif \ --calcA-B \ --NoDataValue-9999 \ --typeFloat32 \ --overwrite为什么要A-B而不是AB回到公式h H N现在已知椭球高h和N求正高H所以H h - N。后面我会专门讲这个最容易搞反的地方。--typeFloat32也建议加上否则遇到整型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:43265773 \ -r bilinear \ dem_ellip.tif \ dem_egm96_via_gdalwarp.tif思路是输入是带椭球高的WGS84三维坐标EPSG:4979输出是WGS84二维加EGM96正高的复合坐标系EPSG:43265773。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 rD:\data\dem_ellip.tif geoid_path rD:\geoid\egm96_15.tif output_path rD:\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), dtypenp.float32) reproject( sourcegeoid_src.read(1), destinationgeoid, src_transformgeoid_src.transform, src_crsgeoid_src.crs, dst_transformsrc.transform, dst_crsdst_crs, resamplingResampling.bilinear, ) dem dem.astype(np.float32) dem[dem profile[nodata]] np.nan corrected dem - geoid profile.update(dtypefloat32, nodata-9999, compressdeflate) 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 projlonglat datumWGS84 geoidgrids/path/to/egm96_15.gtx三是如果你在Python里用pyproj装完pyproj后可以执行pyproj.datadir.get_data_dir()看它到底在找哪个目录然后把网格文件丢进去。这个坑比较烦很多人以为是自己命令写错了其实只是环境问题。4.3 Copernicus DEM是EGM2008怎么转成EGM96Copernicus DEM的垂直基准在公开文档里写得很清楚是EGM2008。如果你非要把EGM2008的高程数据转到EGM96上可以用一个间接公式H96 H2008 N2008 - N96也就是说需要分别算出该点EGM2008的大地水准面差距和EGM96的大地水准面差距然后对高程做一个微调。ICGEM支持同时输出不同模型的N值很方便。操作时先在ICGEM下载EGM2008和EGM96两张网格分别重采样到DEM的范围然后用栅格计算器--calcA B - C其中A是原DEMB是EGM2008网格C是EGM96网格。注意EGM96和EGM2008在全球大部分地区差异在1米以内但部分高山区域可能到3到5米。如果你的项目精度要求不高可能可以忽略但如果涉及水库库容、建筑标高这类精细计算最好转一下。4.4 重采样后边缘出现大片NoData这个问题几乎每次裁剪网格都会遇到。原因很简单原始EGM96网格的范围是全球的而你的DEM范围是局部的gdalwarp在裁剪时由于边界对齐问题可能导致边缘有一圈像元落到了原始网格有效范围外。解决办法有两个。第一个是在gdalwarp里加-wo SAMPLE_STEPS1或者加大-te范围让裁剪框稍微比DEM范围大一圈。第二个更简单做完减法后用原DEM的NoData掩膜把边缘清一遍gdal_calc.py -A dem_egm96.tif -B dem_原始.tif --calcwhere(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值也比收藏夹里吃灰强。
搞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这类文件。虽然这些年EGM20082159阶、EGM20202190阶都出来了EGM96在行业里依然无处不在。原因很简单大量存量数据是以它为垂直基准的。SRTM、NASADEM、ASTER GDEM、ALOS AW3D30这些全球公开DEM产品文档里十有八九都写着“reference to EGM96 geoid”。你如果手里拿着这些数据去做校正反而会把本来正确的值改错。EGM96地位更像是一个“历史基准”新老数据对它都有依赖做数据融合时你绕不开它。1.3 哪些DEM需要校正哪些不需要判断核心就一条看产品元数据里“垂直基准”写的是什么。我按常见数据源整理了一个速查表数据源常见垂直基准是否需要做EGM96校正SRTM1弧秒/3弧秒EGM96不需要直接使用NASADEMEGM96不需要ASTER GDEM v3EGM96不需要ALOS AW3D30EGM96不需要Copernicus DEMGLO-30/GLO-90EGM2008要转EGM96时需转换原始LiDAR点云RTK/PPK解算常见为WGS84椭球高需要校正InSAR处理后的高程产品常见为椭球高需要校正我见过不少朋友拿到Copernicus DEM直接拿来做测区背景没注意它是EGM2008基准结果和其他EGM96数据一并使用后局部区域出现一两米的系统性偏差。所以工作流第一步永远是确认你的数据基准而不是急着跑算法。1.4 工具选型在线查询、网格下载与本地处理校正过程中你一定会用到“查N值”和“批量改栅格”两个操作。查N值适合用在线工具简单快批量校正适合用GDAL或Python脚本稳定可复现。我常用的工具如下工具类型用途入口/查找方式ICGEM Calculation Service在线全球任意点EGM96大地水准面差距计算、网格下载icgem.gfz-potsdam.deNOAA 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:4979WGS84三维地理坐标高程轴默认就是椭球高而且坐标本身是用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_xyTrue) 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 \ --outfiledem_egm96.tif \ --calcA-B \ --NoDataValue-9999 \ --typeFloat32 \ --overwrite为什么要A-B而不是AB回到公式h H N现在已知椭球高h和N求正高H所以H h - N。后面我会专门讲这个最容易搞反的地方。--typeFloat32也建议加上否则遇到整型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:43265773 \ -r bilinear \ dem_ellip.tif \ dem_egm96_via_gdalwarp.tif思路是输入是带椭球高的WGS84三维坐标EPSG:4979输出是WGS84二维加EGM96正高的复合坐标系EPSG:43265773。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 rD:\data\dem_ellip.tif geoid_path rD:\geoid\egm96_15.tif output_path rD:\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), dtypenp.float32) reproject( sourcegeoid_src.read(1), destinationgeoid, src_transformgeoid_src.transform, src_crsgeoid_src.crs, dst_transformsrc.transform, dst_crsdst_crs, resamplingResampling.bilinear, ) dem dem.astype(np.float32) dem[dem profile[nodata]] np.nan corrected dem - geoid profile.update(dtypefloat32, nodata-9999, compressdeflate) 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 projlonglat datumWGS84 geoidgrids/path/to/egm96_15.gtx三是如果你在Python里用pyproj装完pyproj后可以执行pyproj.datadir.get_data_dir()看它到底在找哪个目录然后把网格文件丢进去。这个坑比较烦很多人以为是自己命令写错了其实只是环境问题。4.3 Copernicus DEM是EGM2008怎么转成EGM96Copernicus DEM的垂直基准在公开文档里写得很清楚是EGM2008。如果你非要把EGM2008的高程数据转到EGM96上可以用一个间接公式H96 H2008 N2008 - N96也就是说需要分别算出该点EGM2008的大地水准面差距和EGM96的大地水准面差距然后对高程做一个微调。ICGEM支持同时输出不同模型的N值很方便。操作时先在ICGEM下载EGM2008和EGM96两张网格分别重采样到DEM的范围然后用栅格计算器--calcA B - C其中A是原DEMB是EGM2008网格C是EGM96网格。注意EGM96和EGM2008在全球大部分地区差异在1米以内但部分高山区域可能到3到5米。如果你的项目精度要求不高可能可以忽略但如果涉及水库库容、建筑标高这类精细计算最好转一下。4.4 重采样后边缘出现大片NoData这个问题几乎每次裁剪网格都会遇到。原因很简单原始EGM96网格的范围是全球的而你的DEM范围是局部的gdalwarp在裁剪时由于边界对齐问题可能导致边缘有一圈像元落到了原始网格有效范围外。解决办法有两个。第一个是在gdalwarp里加-wo SAMPLE_STEPS1或者加大-te范围让裁剪框稍微比DEM范围大一圈。第二个更简单做完减法后用原DEM的NoData掩膜把边缘清一遍gdal_calc.py -A dem_egm96.tif -B dem_原始.tif --calcwhere(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值也比收藏夹里吃灰强。
