1. 为什么我要自己动手做ASTER GDEM镶嵌第一次接触ASTER GDEM是在做一个区域地形分析的小项目当时需要覆盖三个省份的连续高程数据。我一开始想得很简单直接去下载现成的分幅数据用GIS软件拼一下不就完了。结果打开一看每幅数据大概1度乘1度的范围光目标区域就涉及几十幅而且相邻图幅之间还有重叠带和边缘噪声。手动拼了几次接边处要么出现明显的亮度断层要么高程值在重叠区对不上做坡度分析的时候直接出现条带状异常。后来我才认真去研究ASTER GDEM这个数据源本身。它是由搭载在Terra卫星上的ASTER传感器获取的立体像对生成的全球数字高程模型第一版在2009年发布第二版在2011年更新空间分辨率是1弧秒也就是大约30米。这个数据覆盖了北纬83度到南纬83度之间的所有陆地对于没有实测地形数据的区域来说几乎是唯一能免费拿到的中高分辨率高程来源。但它的原始分发方式就是按1度乘1度的分幅打包全球一共两万两千多幅你要用哪块就得自己镶嵌。镶嵌这件事听起来不复杂但真正做起来坑不少。首先是数据量的问题一幅GeoTIFF大概二十多兆几十幅加起来就是几个G如果涉及上百幅普通笔记本的内存和磁盘读写就会吃紧。其次是接边处理ASTER GDEM不同图幅之间由于成像时间、大气条件、传感器姿态的差异边缘区域的高程值会有系统性偏差直接拼接会出现明显的接缝。再就是无效值处理水体、云覆盖、陡峭阴影区域会出现-9999这样的填充值如果不做掩膜和插补后续分析全废。我踩过的最大的坑是第一次镶嵌完没有做投影统一。ASTER GDEM原始数据是地理坐标系经纬度以WGS84为基准但如果你拿到的分幅数据里混入了不同投影版本或者你自己在镶嵌过程中不小心改了投影参数拼出来的结果就会出现错位。我当时有一幅数据不知道什么时候被转成了UTM结果那一块整体偏移了几百米做等高线的时候直接穿帮。所以这篇文章我想把整个ASTER GDEM镶嵌的流程从头到尾讲清楚包括数据准备、投影统一、接边处理、无效值填充、镶嵌后裁剪和质检。适合谁看呢我觉得三类人比较合适一是做区域地形分析的学生和研究人员二是做遥感影像处理的技术人员三是需要高程数据做三维建模或者水文分析的朋友。哪怕你之前没碰过GDEM跟着走一遍也能把流程跑通。2. 镶嵌前的整体设计与核心思路拆解2.1 为什么不能直接拼接而要设计流程很多人觉得镶嵌就是把几幅图放到一起软件里点一下拼接就完事了。但ASTER GDEM的分幅数据有几个特殊性质决定了你不能这么干。第一分幅之间是有重叠的。ASTER GDEM的每一幅标准分幅是1度乘1度但实际有效数据范围会略大于这个边界相邻图幅之间通常有几十个像素的重叠带。这个重叠带的存在本来是为了方便接边但如果你直接拼接重叠区域就会出现两套高程值软件默认取其中一个接缝处就会出现突变。第二不同图幅的高程基准有微小差异。虽然ASTER GDEM整体是基于WGS84椭球的高程但不同图幅在生成过程中使用的控制数据不同导致图幅之间可能存在米级的系统性偏差。这个偏差在平坦地区不明显但在坡度分析、水文流向计算中会被放大。第三无效值分布不均匀。水体、云、雪、阴影区域都会产生无效值而且不同图幅的无效值分布不一样。如果你先镶嵌再处理无效值会发现有些接缝处的无效值连成一片很难判断是数据本身的问题还是拼接引入的。所以我的整体思路是先做单幅数据的预处理包括投影统一、无效值标记、边缘噪声裁剪然后做接边处理通过重叠区的高程统计来调整图幅之间的偏差最后才做镶嵌和裁剪。这个顺序不能乱乱了就要返工。2.2 工具选型为什么我最终用了GDAL加Python工具选型这块我试过不少方案。ArcGIS的镶嵌工具确实方便Mosaic To New Raster点几下就能出结果但问题是它对大批量数据的处理效率不高而且接边处理的参数不够透明你很难控制它到底怎么取重叠区的值。QGIS的虚拟栅格功能适合做快速预览但真正要输出一个干净的结果还是得靠命令行工具。我最终选的是GDAL加Python脚本的方案。GDAL的gdalwarp可以做投影转换和重采样gdalbuildvrt可以快速构建虚拟镶嵌gdal_translate可以把虚拟镶嵌转成实体文件。Python的rasterio和numpy可以让我灵活地读取重叠区的高程值做统计分析和偏差调整。这个组合的好处是每一步都可控参数都在你手里出了问题也能定位到具体是哪一步。提示如果你只是做一次性的小区域镶嵌用QGIS的虚拟栅格加导出就够了。但如果涉及几十幅以上的数据或者需要反复调整接边参数建议还是走脚本路线。2.3 投影选择与重采样方法ASTER GDEM原始数据是地理坐标系单位是度。做镶嵌的时候我建议统一转成投影坐标系单位是米。原因很简单地理坐标系下不同纬度的1度对应的实际距离不一样你在做坡度、距离、面积计算的时候会引入纬度相关的误差。转成投影之后单位统一计算也方便。投影的选择要看你的研究区域。如果区域不大比如一个省的范围用UTM分区投影就够了。如果区域跨多个UTM带可以用Albers等面积投影或者Lambert等角投影。我一般用Albers因为它在中等纬度地区的面积变形小适合做区域分析。重采样方法我试过最近邻、双线性、三次卷积。最近邻速度快但会引入阶梯状的高程突变做坡度分析的时候会出现明显的块状噪声。双线性平滑一些但在陡峭地形区域会过度平滑。三次卷积效果最好但计算量大。我的建议是如果做坡度、水文分析用双线性如果做三维可视化用三次卷积如果只是做粗略预览最近邻也行。# 投影转换示例将地理坐标系转为Albers投影 gdalwarp -t_srs projaea lat_125 lat_247 lat_036 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 30 30 -r bilinear \ input_gdem.tif output_gdem_proj.tif这段命令里-t_srs指定目标投影-tr 30 30指定输出分辨率为30米-r bilinear指定双线性重采样。注意-tr后面的两个30是x和y方向的分辨率单位是米和投影单位一致。3. 核心细节解析与实操要点3.1 数据下载与命名规范ASTER GDEM的数据可以从多个渠道获取我常用的是NASA的EarthData平台。下载的时候要注意每一幅数据都有唯一的标识码命名格式一般是ASTGTMV003_NxxExxx这样的结构其中N后面跟纬度E后面跟经度。下载下来是压缩包解压后有两个文件一个是dem文件一个是num文件。dem就是高程数据num是质量评估文件记录了每个像素的有效性。我建议在下载阶段就做好命名规范。比如你研究区域是北纬30到33度东经100到103度那涉及的分幅就是N30E100到N32E102一共九幅。你可以按行列号重命名比如row30_col100.tif这样后面做镶嵌的时候排序和索引都方便。注意下载的时候一定要确认版本号。ASTER GDEM有V1、V2、V3三个版本V3是最新的修正了V2中的一些异常值。如果你拿到的数据里混了不同版本接边处的高程偏差会更大。3.2 无效值处理与质量掩膜ASTER GDEM的无效值通常用-9999表示。这些无效值出现在水体、云覆盖、传感器阴影、陡峭地形等区域。如果你不做处理镶嵌后的结果里就会有大片的-9999做任何分析都会出错。我的处理方法是分两步。第一步用num文件生成质量掩膜。num文件里每个像素的值代表了该像素的有效性等级一般大于0的值表示有效等于0表示无效。你可以用gdal_calc把num文件转成二值掩膜然后应用到dem上。第二步对于掩膜后的无效区域用邻域插值或者全局插值填充。邻域插值适合小面积的无效区全局插值适合大面积的缺失。import rasterio import numpy as np from scipy import ndimage # 读取dem和num文件 with rasterio.open(dem.tif) as src: dem src.read(1) profile src.profile with rasterio.open(num.tif) as src: num src.read(1) # 生成有效掩膜 valid_mask (num 0) (dem ! -9999) # 对无效区域做邻域均值填充 filled_dem dem.copy() filled_dem[~valid_mask] np.nan # 用3x3窗口的均值填充 kernel np.ones((3,3)) / 9 neighbor_mean ndimage.convolve(np.nan_to_num(filled_dem), kernel, modeconstant) filled_dem[~valid_mask] neighbor_mean[~valid_mask] # 如果还有剩余的nan用全局均值填充 global_mean np.nanmean(filled_dem) filled_dem[np.isnan(filled_dem)] global_mean # 写出结果 with rasterio.open(dem_filled.tif, w, **profile) as dst: dst.write(filled_dem, 1)这段代码的逻辑是先用num文件生成有效掩膜然后把无效区域设为nan用3x3窗口的邻域均值填充最后用全局均值兜底。实际用的时候窗口大小可以根据无效区的面积调整大面积缺失可以用更大的窗口或者用反距离加权插值。3.3 接边处理的核心逻辑接边处理是镶嵌里最关键的环节。我的做法是先找到相邻图幅的重叠区然后在重叠区里统计两幅图的高程差异如果差异超过阈值就对其中一幅做整体偏移调整。具体步骤是这样的。假设你有两幅图A和B它们在水平方向上有重叠。你先用gdalwarp把两幅图都转到同一个投影和分辨率下然后用gdalbuildvrt构建一个虚拟镶嵌再用Python读取重叠区的像素值。计算A和B在重叠区的平均高程差如果差值大于5米就把B的整体高程加上这个差值。这个5米的阈值是我根据经验定的平坦地区可以设小一点山区可以设大一点。提示接边调整不要一次性对所有图幅做建议按行列顺序逐对处理。先处理第一行再处理第二行每一行内部从左到右依次调整。这样可以把累积误差控制住。3.4 边缘噪声裁剪ASTER GDEM的每一幅数据边缘都有大约1到2个像素的噪声带表现为高程值的突变或者无效值。如果你不裁剪镶嵌后的接缝处就会出现细条状的异常。我的做法是在镶嵌前对每一幅数据做边缘裁剪裁掉最外圈的2个像素。# 边缘裁剪示例裁掉最外圈2个像素 gdal_translate -srcwin 2 2 $(expr $(gdalinfo input.tif | grep Size is | awk {print $3}) - 4) $(expr $(gdalinfo input.tif | grep Size is | awk {print $4}) - 4) input.tif output_trimmed.tif这个命令用gdal_translate的-srcwin参数指定裁剪窗口从第2行第2列开始宽度和高度各减去4个像素。实际用的时候你可以先打开一幅数据看看边缘噪声的宽度再决定裁多少。4. 完整实操流程与关键环节实现4.1 环境准备与依赖安装我用的环境是Ubuntu 20.04加Python 3.8。GDAL的安装建议用conda因为pip安装的GDAL经常和系统库版本不匹配。rasterio和numpy也用conda装scipy用来做插值。conda create -n gdem python3.8 conda activate gdem conda install -c conda-forge gdal rasterio numpy scipy装完之后验证一下GDAL版本确保gdalwarp和gdalbuildvrt都能正常调用。gdalinfo --version如果输出类似GDAL 3.4.1这样的版本信息就说明装好了。4.2 批量投影转换与重采样假设你已经下载了九幅数据放在data目录下命名是N30E100.tif到N32E102.tif。第一步是批量做投影转换。#!/bin/bash for file in data/*.tif; do basename$(basename $file .tif) gdalwarp -t_srs projaea lat_125 lat_247 lat_036 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 30 30 -r bilinear \ $file proj/${basename}_proj.tif done这个脚本会遍历data目录下的所有tif文件逐个做投影转换输出到proj目录。注意-tr 30 30指定了输出分辨率这个值要和你的研究需求匹配。如果你需要更高分辨率可以设成15或者10但那样数据量会翻倍。4.3 无效值填充与边缘裁剪投影转换完之后对每一幅数据做无效值填充和边缘裁剪。我把这两步合并到一个Python脚本里。import rasterio import numpy as np from scipy import ndimage import os def fill_and_trim(input_path, output_path, trim_pixels2): with rasterio.open(input_path) as src: dem src.read(1) profile src.profile # 无效值填充 valid_mask dem ! -9999 filled_dem dem.copy() filled_dem[~valid_mask] np.nan # 邻域均值填充 kernel np.ones((3,3)) / 9 neighbor_mean ndimage.convolve(np.nan_to_num(filled_dem), kernel, modeconstant) filled_dem[~valid_mask] neighbor_mean[~valid_mask] # 全局均值兜底 global_mean np.nanmean(filled_dem) filled_dem[np.isnan(filled_dem)] global_mean # 边缘裁剪 trimmed_dem filled_dem[trim_pixels:-trim_pixels, trim_pixels:-trim_pixels] # 更新profile profile.update({ height: trimmed_dem.shape[0], width: trimmed_dem.shape[1], transform: rasterio.transform.from_origin( src.transform.c trim_pixels * src.transform.a, src.transform.f trim_pixels * src.transform.e, src.transform.a, -src.transform.e ) }) with rasterio.open(output_path, w, **profile) as dst: dst.write(trimmed_dem, 1) # 批量处理 for file in os.listdir(proj): if file.endswith(_proj.tif): input_path os.path.join(proj, file) output_path os.path.join(filled, file.replace(_proj.tif, _filled.tif)) fill_and_trim(input_path, output_path)这个脚本里trim_pixels2表示裁掉最外圈2个像素。rasterio.transform.from_origin用来更新裁剪后的地理变换参数确保裁剪后的数据还能正确定位。4.4 接边偏差统计与调整接边调整我写了一个单独的脚本逻辑是遍历所有相邻图幅对计算重叠区的高程差然后做偏移调整。import rasterio import numpy as np import os def get_overlap_diff(file_a, file_b): with rasterio.open(file_a) as src_a: dem_a src_a.read(1) transform_a src_a.transform crs_a src_a.crs with rasterio.open(file_b) as src_b: dem_b src_b.read(1) transform_b src_b.transform crs_b src_b.crs # 确保投影一致 if crs_a ! crs_b: return None # 计算重叠区 # 这里简化处理假设两幅图在水平方向相邻 # 实际用的时候需要用rasterio的window功能精确计算 overlap_width min(dem_a.shape[1], dem_b.shape[1]) // 10 overlap_a dem_a[:, -overlap_width:] overlap_b dem_b[:, :overlap_width] # 计算平均高程差 diff np.nanmean(overlap_a) - np.nanmean(overlap_b) return diff # 遍历相邻图幅 files sorted([f for f in os.listdir(filled) if f.endswith(.tif)]) for i in range(len(files) - 1): diff get_overlap_diff( os.path.join(filled, files[i]), os.path.join(filled, files[i1]) ) if diff is not None and abs(diff) 5: print(f调整 {files[i1]}偏移量{diff:.2f}米) # 实际调整逻辑读取文件加上偏移量重新写出这个脚本是个简化版实际用的时候需要用rasterio的window功能精确计算重叠区而不是简单取边缘。但核心逻辑是一样的找到重叠区算平均差超过阈值就调整。4.5 虚拟镶嵌与实体输出所有单幅数据处理好之后用gdalbuildvrt构建虚拟镶嵌再用gdal_translate转成实体文件。# 构建虚拟镶嵌 gdalbuildvrt mosaic.vrt filled/*.tif # 转成实体GeoTIFF gdal_translate -of GTiff -co COMPRESSLZW -co TILEDYES mosaic.vrt mosaic.tif-co COMPRESSLZW指定用LZW压缩可以显著减小文件体积。-co TILEDYES指定用分块存储方便后续读取。4.6 镶嵌后裁剪与质检镶嵌完的mosaic.tif覆盖了你下载的所有图幅范围但你可能只需要其中一部分。用gdalwarp加-cutline或者用gdal_translate加-projwin做裁剪。# 按经纬度范围裁剪 gdal_translate -projwin 100.0 33.0 103.0 30.0 mosaic.tif final.tif-projwin后面跟的是左上角经度、左上角纬度、右下角经度、右下角纬度。注意这个顺序不能错否则裁出来是反的。质检我一般做三件事。第一用gdalinfo看统计信息确认最大值、最小值、均值在合理范围内。第二用QGIS打开看接缝处有没有明显的条带。第三随机取几个点和已知的高程数据做对比。提示质检的时候特别要注意接缝处的高程突变。如果发现某条接缝有明显的亮度差异说明接边调整没做好需要回去重新处理那一对图幅。5. 常见问题与排查技巧实录5.1 接缝处出现条带状异常怎么办这是最常见的问题。原因通常有三个一是接边偏差没调整二是边缘噪声没裁剪三是重采样方法不一致。排查顺序是这样的。先检查所有图幅的投影和分辨率是否一致用gdalinfo逐个看。如果投影不一致回去重新做投影转换。如果投影一致检查边缘裁剪是否都做了裁了多少像素。如果都没问题那就是接边偏差的问题需要重新计算重叠区的高程差调整阈值。我遇到过一次接缝处出现了一条宽度约5个像素的亮带。查了半天发现是其中一幅数据在做投影转换的时候用了最近邻重采样其他用的是双线性。最近邻在边缘产生了阶梯状的高程突变和双线性的平滑结果拼在一起就出现了亮带。统一重采样方法之后就解决了。5.2 镶嵌后数据量太大怎么优化几十幅数据镶嵌完文件可能有好几个G。优化方法有几个。第一用LZW或者DEFLATE压缩一般能压到原来的三分之一。第二用分块存储TILEDYES这样读取的时候不用加载整个文件。第三如果研究区域不大先裁剪再镶嵌而不是先镶嵌再裁剪。第四如果不需要原始分辨率可以重采样到更低的分辨率比如从30米降到90米数据量直接降到九分之一。# 压缩和分块存储 gdal_translate -of GTiff -co COMPRESSDEFLATE -co PREDICTOR2 -co TILEDYES -co BIGTIFFIF_NEEDED mosaic.tif mosaic_compressed.tif-co PREDICTOR2对高程数据特别有效因为它利用了相邻像素的相关性。-co BIGTIFFIF_NEEDED确保文件超过4G的时候自动用BigTIFF格式。5.3 无效值填充后出现异常高值或低值这个问题通常出现在大面积无效区填充的时候。邻域均值填充在边缘区域会引入异常值因为邻域里可能包含无效值或者极端值。我的解决方法是先用有效掩膜把无效区标记出来然后只用有效像素做邻域统计。如果无效区面积超过窗口大小就扩大窗口或者改用反距离加权插值。还有一个技巧是用中值滤波代替均值滤波。中值滤波对极端值不敏感填充结果更稳健。但中值滤波的计算量比均值滤波大大面积填充的时候会慢一些。5.4 常见问题速查表问题现象可能原因排查方法解决方案接缝处条带状异常接边偏差未调整检查重叠区高程差重新计算偏差并调整接缝处细线状异常边缘噪声未裁剪检查边缘像素值裁剪最外圈2-3像素镶嵌后文件过大未压缩或未分块查看文件属性用LZW压缩和TILED存储无效值填充后异常邻域包含无效值检查填充区域用有效掩膜限制邻域投影错位投影未统一gdalinfo检查CRS统一投影后重新镶嵌坡度分析出现块状噪声重采样方法不当检查重采样参数改用双线性或三次卷积5.5 几个我踩过的坑和独家技巧第一个坑是下载数据的时候没注意版本号。我有一次混用了V2和V3的数据接边处的高程差达到了十几米调了半天才发现是版本问题。后来我养成了习惯下载的时候先看文件名里的版本标识统一用V3。第二个坑是投影转换的时候没设分辨率。gdalwarp默认会保持原始分辨率但地理坐标系下的1度对应的实际距离和投影坐标系下的30米不是一回事。如果你不设-tr参数转换后的分辨率可能是乱的。我一般显式设成30 30确保输出分辨率一致。第三个技巧是接边调整的时候用中位数而不是平均值。平均值容易受极端值影响中位数更稳健。特别是在山区重叠区里可能有陡峭地形导致的高程突变用中位数可以避免这些突变影响整体偏移量。第四个技巧是镶嵌前先做一个低分辨率的预览。用gdalbuildvrt加-resolution lowest参数快速生成一个低分辨率镶嵌看看整体效果。如果预览里就有明显的接缝那高分辨率的结果肯定也有问题可以提前发现。第五个技巧是保留中间文件。投影转换后的文件、填充后的文件、裁剪后的文件都保留着。万一后面发现某一步有问题不用从头再来。磁盘空间不够的话至少保留投影转换后的文件因为那一步最耗时。6. 镶嵌后的数据应用与扩展思路6.1 基于镶嵌结果的地形分析镶嵌完的GDEM数据可以直接用来做坡度、坡向、地形起伏度分析。我用GDAL的gdaldem工具做坡度分析一条命令就能出结果。# 坡度分析 gdaldem slope final.tif slope.tif -p -s 1.0 # 坡向分析 gdaldem aspect final.tif aspect.tif # 地形起伏度用3x3窗口的极差 gdaldem roughness final.tif roughness.tif-p表示输出以度为单位的坡度-s 1.0表示垂直方向的缩放系数。如果你的高程数据单位是米投影单位也是米这个系数就是1.0。如果投影单位是度就需要调整这个系数。坡度分析的结果可以用来做地形分级、水土流失评估、选址分析等。我做过一个区域的地形起伏度分析用镶嵌后的GDEM算出来的结果和实地考察的数据对比误差在可接受范围内。6.2 水文分析中的流向计算GDEM镶嵌结果在水文分析里用得很多比如计算流向、汇流累积量、河网提取。但GDEM本身有噪声直接算流向会出现很多虚假的洼地。我的做法是先做洼地填充用GDAL的fillnodata或者用Python的richdem库。import richdem as rd # 读取镶嵌后的DEM dem rd.LoadGDAL(final.tif) # 洼地填充 dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) # 计算流向 flow_dir rd.FlowDirD8(dem_filled, in_placeFalse) # 计算汇流累积量 flow_acc rd.FlowAccumulation(flow_dir, methodD8) # 保存结果 rd.SaveGDAL(flow_acc.tif, flow_acc)richdem的FillDepressions函数可以自动填充洼地epsilonTrue表示用微小坡度填充避免平坦区域出现流向不确定的问题。汇流累积量可以用来提取河网设定一个阈值比如1000个像素超过这个阈值的区域就是河道。6.3 三维可视化与场景搭建镶嵌后的GDEM可以导入Blender或者Unity做三维地形可视化。我一般把GDEM转成PNG的高度图然后在Blender里用Displace修改器生成地形。具体做法是先用gdal_translate把GDEM转成16位PNG然后在Blender里创建一个平面细分足够多的网格用Displace修改器加载PNG作为高度图。提示转PNG的时候要注意拉伸对比度。GDEM的高程范围可能很大直接转PNG会丢失细节。我一般先用gdal_calc做线性拉伸把高程范围映射到0到65535之间再转PNG。6.4 后续可以扩展的方向这个镶嵌流程还可以扩展到其他GDEM数据源比如SRTM、ALOS World 3D。不同数据源的镶嵌逻辑类似但接边处理的参数需要调整。SRTM的分辨率是90米ALOS World 3D是30米精度和覆盖范围也不一样。如果你需要更高精度的数据可以考虑用商业数据源但那个成本就上去了。另一个扩展方向是自动化。我现在是把整个流程拆成几个脚本手动按顺序跑。如果你需要频繁做镶嵌可以把这些脚本串起来加一个配置文件指定输入目录、输出目录、投影参数、分辨率、裁剪范围一键跑完。我用Python的argparse写了一个简单的命令行工具支持这些参数用起来方便不少。还有一个方向是质量评估的自动化。我现在质检是靠人工看效率不高。后来我写了一个脚本自动计算镶嵌后数据的统计信息包括均值、标准差、无效值比例、接缝处的高程梯度。如果某个指标超过阈值就报警提示。这个脚本帮我省了不少时间特别是做大批量镶嵌的时候。我个人在实际操作中的体会是ASTER GDEM镶嵌这件事难点不在技术本身而在细节的把控。投影要统一分辨率要一致接边要调整无效值要填充边缘要裁剪每一步都不能省。省了一步后面就要花更多时间返工。另外中间文件一定要保留特别是投影转换后的文件因为那一步最耗时万一后面发现问题不用从头再来。最后再分享一个小技巧如果你不确定接边调整的阈值设多少合适可以先设一个保守的值比如3米跑一遍看看效果如果接缝还是明显再逐步放宽到5米、10米。不要一上来就设很大的阈值那样可能会把真实的地形差异也抹掉。