4. DEM预处理：数据格式转换、投影变换、空洞填充与滤波
各位同学，咱们今天聊DEM预处理。说实话，很多人拿到DEM数据就直接往软件里一扔，开始做三维建模。我当年也这么干过，结果模型出来全是锯齿，地形分析根本没法看。后来才明白——预处理这一步，决定了你后面所有工作的质量。
DEM预处理说白了就四件事：格式转换、投影变换、空洞填充、滤波去噪。咱们一个一个来拆解。
4.1 数据格式转换——别让格式卡住你
DEM数据的格式五花八门。我见过最头疼的是有人拿了个.asc文件，愣是不知道怎么打开。其实常见的格式就那么几种：
| 格式 | 特点 | 适用场景 |
|---|---|---|
| GeoTIFF | 自带地理信息，通用性强 | 绝大多数GIS软件 |
| ASCII Grid (.asc) | 纯文本，可读性好 | 小范围数据交换 |
| IMG | ERDAS格式，支持多波段 | 遥感专业软件 |
| DTED | 军事标准格式 | 国防、航空领域 |
我个人习惯，所有DEM数据统一转成GeoTIFF。为什么？因为它把坐标信息、投影信息都嵌在文件里，不会丢。我曾经接手过一个项目，对方给了几十个.asc文件，结果投影信息全写在单独的文本里，对不上号，折腾了两天才对齐。
核心原则：格式转换不是简单的改后缀名。一定要确保转换过程中不丢失高程精度，不改变数据范围。
用GDAL做格式转换，一行命令搞定：
gdal_translate -of GTiff input.asc output.tif嗯，就这么简单。但要注意，如果源数据是浮点型，转成整型会损失精度。我建议始终保留浮点型输出。
4.2 投影变换——坐标系不对，一切白费
投影变换是DEM预处理里最容易出问题的地方。你想想看，两个不同投影的数据叠在一起，位置差了几百米，那地形分析还有什么意义？
我遇到过最离谱的事：有人把WGS84经纬度的DEM直接当成平面坐标用，算出来的坡度全是错的。后来查了半天才发现是投影没转。
常见的投影坐标系：
- UTM投影——适合中小范围，变形小
- 高斯-克吕格投影——国内常用，按3度带或6度带分
- Albers等面积投影——适合面积计算
我的建议：如果做地形分析，优先用UTM投影。因为它保持角度和形状，坡度、坡向计算最准确。
用GDAL做投影变换：
gdalwarp -t_srs EPSG:32650 input.tif output_utm.tifEPSG:32650是UTM 50N的编码。你查一下你所在区域的UTM带号，别搞错了。我曾经帮一个学生调试，他用了32649，结果数据偏了半个中国。
注意：投影变换时，重采样方法会影响高程值。我一般用三次卷积（cubic），精度高。如果追求速度，可以用双线性（bilinear）。千万别用最近邻（nearest），那会让地形出现台阶状。
4.3 空洞填充——数据缺失怎么办？
DEM数据经常有空洞。原因很多：云遮挡、雷达阴影、传感器故障。空洞不补，后面做水文分析、可视域分析都会出问题。
空洞填充的方法，我按推荐程度排个序：
- 插值法——用周围像素推算空洞值。最常用的是反距离权重法（IDW）和克里金法。
- 填补法——用其他数据源填补。比如SRTM的空洞用ASTER GDEM补。
- 平滑法——对空洞边缘做平滑过渡，适合小空洞。
我个人习惯，小空洞（小于10个像素）用插值法，大空洞用填补法。为什么？因为插值法在大空洞区域会产生不自然的平坦区域，你想想看，一个山谷突然变成平地，那地形分析就失真了。
用ArcGIS做空洞填充：
# 使用焦点统计工具
Focal Statistics (Con, 3x3矩形窗口, 平均值)或者用GDAL的fillnodata工具：
gdal_fillnodata.py input.tif output_filled.tif -md 10-md 10表示最大搜索距离为10个像素。这个值要根据空洞大小调整。我曾经处理过一个LiDAR数据，空洞直径有50个像素，md设成10根本填不上，后来改成30才搞定。
避坑指南：空洞填充后一定要检查。我曾经填完没检查，结果在河流区域填出了一个小山包，水文分析全错了。后来我养成了习惯——填充后做一次坡度分析，看看有没有异常突变。
4.4 滤波去噪——让地形更平滑
DEM数据里总有噪声。可能是传感器误差，也可能是地物干扰（比如树木、建筑物）。滤波的目的就是去掉这些噪声，保留真实地形。
常用的滤波方法：
| 滤波方法 | 原理 | 适用场景 |
|---|---|---|
| 均值滤波 | 取窗口内平均值 | 去除随机噪声 |
| 中值滤波 | 取窗口内中位数 | 去除椒盐噪声，保留边缘 |
| 高斯滤波 | 加权平均，中心权重高 | 平滑效果好，但会模糊边缘 |
| 双边滤波 | 同时考虑距离和灰度差异 | 保留边缘的同时去噪 |
嗯，这里要注意。均值滤波虽然简单，但会把地形细节也抹掉。我做过对比实验，3x3的均值滤波，山脊线会变模糊，坡度值平均降低5%-10%。
我个人推荐中值滤波。为什么？因为它能有效去除孤立噪声点，同时保留地形边缘。你想想看，一个突然的尖峰（可能是树冠），中值滤波直接把它当成异常值去掉，而均值滤波会把它扩散到周围。
用GDAL做中值滤波：
gdal_filter.py input.tif output_filtered.tif -filter median -radius 3半径3表示3x3窗口。如果噪声严重，可以试试5x5，但要注意地形细节也会损失更多。
我的经验：滤波窗口大小要根据DEM分辨率来定。30米分辨率的DEM，3x3窗口覆盖90米范围，足够了。如果是5米分辨率的高精度DEM，可以试试5x5甚至7x7。
4.5 预处理流程总结
说了这么多，咱们把整个流程串起来。我画了一张图，你看一眼就明白了：
你看，整个流程是线性的。但要注意，每一步都可能需要回头调整。比如滤波后发现空洞没填好，就得回去重新填。我做过一个项目，来回迭代了三次才满意。
最后说一句：预处理不是走过场。你花在预处理上的每一分钟，都会在后面的地形分析中加倍回报。我见过太多人跳过预处理直接建模，结果模型精度差、分析结果不可靠，最后还得回头重做。嗯，别走弯路。
公众号:蓝海资料掘金营，微信deep3321