今日分享：

DeltaDTM 全球海岸数字高程模型

DeltaDTM 是一款全球海岸数字地形模型 (DTM)，水平空间分辨率为 1 弧秒（约 30 米），垂直平均绝对误差 (MAE) 为 0.43 米。它利用 ICESat-2 和 GEDI 任务的星载激光雷达数据校正哥白尼数字高程模型 (DEM)，从而提高了现有全球高程数据集的精度。校正过程包括偏差校正、非地形单元（例如植被和建筑物）过滤以及使用插值法填充数据缺失。DeltaDTM 特别关注低洼沿海地区（海拔低于平均海平面 30 米），这些地区极易受到海平面上升、地面沉降和极端天气事件的影响。

DeltaDTM 是一项极具价值的数据，可用于包括海岸管理、洪水建模和适应性规划在内的广泛应用。其更高的精度能够更精确地评估海岸洪水风险，并支持制定有效的缓解和适应策略。

文章网址：

https://www.nature.com/articles/s41597-024-03091-9

海岸高程数据对于海岸管理、洪水建模和适应性规划等众多应用至关重要。低洼海岸地区（海拔低于平均海平面10米）面临未来极端水位、地面沉降和极端天气模式变化的风险。然而，目前可免费获取的高程数据集精度不足以模拟这些风险。研究团队推出DeltaDTM，这是一个全球海岸数字地形模型（DTM），可在公共领域获取，其水平空间分辨率为30m，垂直平均绝对误差（MAE）为0.45米。DeltaDTM利用ICESat-2和GEDI任务的星载激光雷达数据对CopernicusDEM进行校正。校正的具体内容有：1.校正了CopernicusDEM中的高程偏差，应用滤波算法去除非地形单元，并使用插值法填充缺失数据。值得注意的是，其团队提出的分类方法比近期其他研究者用于校正数字高程模型（DEM）的回归方法更为精确，后者最佳的平均绝对误差（MAE）也仅为0.72米。

DeltaDTM 是一个基于 CopernicusDEM、ICESat-2 和 GEDI 高程数据融合的全球海岸数字地形模型 (DTM)。利用 ICESat-2 和 GEDI地形高程测量数据，消除了 CopernicusDEM 中存在的地表数据（例如，冠层、建筑物）的垂直偏差。我们的方法可以分为四类：

空间滤波，例如去除 CopernicusDEM 中存在的坑洞和其他异常值。
配准，使 CopernicusDEM 和 ICESat-2 垂直对齐，从而消除 CopernicusDEM 中的垂直偏差。
^{通过使用形态学滤波器 CopernicusDEM 分类为地形和非地形，并去除非地形高程像素，从而过滤非地面点。}
使用 AIDW对上一步移除的值进行空间插值，从而填充空隙。

图 ( a ) 和 ( b ) 分别展示了加里曼丹和荷兰DeltaDTM 的分类过程。顶行显示了CopernicusDEM（DeltaDTM 的输入 DSM）以及该区域的参考机载激光雷达 DTM。中间行显示了以ESA WorldCover地图为参考的地形像素分类结果。底行显示了DeltaDTM ，即地形插值的结果，并以归一化 DSM为参考。归一化 DSM 是通过从 CopernicusDEM 中减去 DeltaDTM 而生成的，最终得到地表高程图。与其他校正后的数字表面模型（例如FABDEM）的不同之处在于，我们使用分类而非回归来确定地形高度。虽然这在理论上会导致分辨率损失（数据经过滤波和插值后），但我们发现，精度的提高可以弥补这一缺陷，尤其是在数据稀少的区域。

使用的数据集：

我们用作基准高程模型的是 CopernicusDEM GLO-30 数据集，该数据集由欧盟和欧洲航天局 (ESA) 在 COPERNICUS 计划下提供，该数据集以 1 度 × 1 度的图块形式分布，空间分辨率为 1 角秒（赤道附近约为 30 米）。它基于 TanDEM-X 干涉合成孔径雷达 (SAR) 数据，每个高程图块都包含水体掩膜和高程误差图块，我们在分析中也使用了这些信息。

为了获得垂直方向精度更高但分布稀疏的地形高程测量数据，我们使用了ICESat-2 3级陆地和植被高度（ATL08）产品，日期范围为2018年10月14日至2023年6月22日。我们从NSIDC DAAC下载了262807个数据块（总计约22 TB）。对于高程，我们使用了h_te_best_fit_20 m（20米范围内所有地形光子的最佳拟合）字段，该字段包含高于WGS84椭球体的高度以及HDF5文件中每个轨迹组对应的纬度latitude_20m和经度longitude_20m字段。

同样，我们下载了 74815 个全球生态系统动力学调查 (GEDI) 2A 级产品数据^{总计约 107 TB ，日期范围从 2019 年 4 月 18 日到 2023 年 3 月 16 日。我们使用HDF5 文件中每个轨迹组的}elev_lowestmode字段，该字段包含高于 WGS84 椭球体的地形高程以及相关的纬度lat_lowestmode和经度lon_lowestmode字段。

^{我们还从 ESA WorldCover 2021 数据集}中采样了土地覆盖类型，以便进一步用于偏差校正和分类算法。选择该土地覆盖类型数据集是因为它数据较新，且分辨率约为 10 米，高于 CopernicusDEM。该数据以 3°×3° 的图块形式自由分发，我们将其重采样到与 CopernicusDEM 相同的图块规格。WorldCover 识别出多种土地覆盖类型，例如“草地”、“耕地”、“树木覆盖”和“人工建成区”。具体而言，在我们的偏差校正和滤波过程中，我们将“灌木丛”、“草地”、“耕地”、“裸地”、“苔藓”和“雪”定义为开放土地覆盖（即未被木本植被或建筑物覆盖的地形），而将剩余的“树木覆盖”、“红树林”和“建成区”定义为封闭土地覆盖。我们假设 CopernicusDEM 中开阔地表覆盖的高程值可以近似测量地形，而封闭地表覆盖的高程值则不能。

数据预处理

我们通过与 GLL_DTM交集以及对洪泛区高程的人工检查，找到了所有包含低于 10 米（高于平均海平面）值的 CopernicusDEM 图块。最终得到 7146 个图块的子集，总共有 26448 个图块。

我们按照 CopernicusDEM 子集中的图块规范对所有数据集进行图块划分。这简化了处理流程，并实现了图块的并行处理。为了防止图块边界出现边缘伪影，每个图块的处理都与相邻图块保持 5% 的重叠。此重叠百分比是一个安全裕度，因为后续的滤波和空隙填充可能需要 12 公里的缓冲区。

空间滤波

CopernicusDEM、ICESat-2 和 GEDI 数据集都包含异常值。高程高于实际地形的异常值会在后续分类步骤中被自动移除（如同它们被分类为建筑物或树冠），而高程低于实际地形的异常值则会造成问题，因为它们也会被分类为地形，从而对周围像元的分类产生负面影响。

CopernicusDEM 中包含许多低位异常值，这通常是城市地区（例如电线杆周围）多次反射后向散射误差造成的。我们应用一个 25×25 像素（约 750×750 米）的窗口函数，并移除窗口内所有低于所有高程值两个标准差的值。该窗口大小足以过滤掉 CopernicusDEM 中观察到的较大区域（3×3 像素）的低位异常值。此外，对于每个 1×1 像素的图块，我们确定一个高程截止值，低于该值的所有数据都将被移除。默认值设置为 -2 米，对于地势最低的区域（例如荷兰的圩田），我们手动调整了截止值，将其设置为 -7 米。该值作为DEM 附带的tiles.gpkg地理空间数据库中的low_cutoff字段提供。由此产生的空白区域将按照后续步骤中的描述进行填充，并使用重叠图块来防止出现任何边缘伪影。

基于CopernicusDEM 提供的高度误差数据移除了所有高度误差超过 0.75 米或 DEM 被其他 DEM 填充的高程值。这些值是根据验证区域中观察到的异常值经验性地选择的。CopernicusDEM 中约有 1% 的值通过这种方式被移除。我们注意到，FABDEM 13也尝试移除 CopernicusDEM 中的低异常值，但它是通过反复平滑高程值来实现的。

星载激光雷达

对 ICESat-2 和 GEDI 数据都应用了质量过滤。对于 ICESat-2 数据，我们仅保留subset_te_flag标志设置为 1 的数据。对于 GEDI 数据，我们应用了与衍生 GEDI L3A 产品相同的过滤方法，并且仅保留sensitivity标志高于 0.95 的数据。与 CopernicusDEM 使用的过滤方法类似，我们也移除了约 750 × 750 米窗口内所有高程值低于 2 个标准差的所有值。此外，我们还使用了相同的高程截止值（默认设置为 -2 米），低于该值的所有数据均被移除。这些过滤方法会移除约 1% 的测量数据。

^{ICESat-2 和 GEDI 高程值均使用 PROJ 30}从椭球体垂直变换到地球重力模型 (EGM2008) 大地水准面^{（在本研究中假设近似于全球平均海平面），与 CopernicusDEM 使用的垂直参考相同}。

CopernicusDEM——现有的全球雷达或光学数字高程模型（DEM）一样测量地球表面，因此包含了植被、建筑物高度。为了消除这些偏差并确定真实的“裸地”表面，我们应用了形态学表面滤波器，这些滤波器由ICESat-2 ATL08和GEDI L2A数据的地形测量数据支持。形态学滤波器与地物形态（形状）相关，并作用于栅格（图像）数据的子区域（窗口），在这些子区域上应用非线性滤波器（例如最小值滤波器）。这些滤波器常用于机载激光雷达数据集的地形分类。

利用ICESat-2 ATL08和GEDI L2A地形数据替换CopernicusDEM数据，将激光雷达获取的高程值“替换”到经过偏差校正的CopernicusDEM栅格数据中。

空隙填充

反距离加权 (AIDW) 方法，利用周围的地形点进行插值填充。该方法是一种标准的反距离加权 (IDW) 方法，但对于相对于插值点“靠后”的点，其权重较低。实际上，这确保了所使用的值来自目标点周围的所有点，而不仅仅是同一方向上的邻近点。这一特性避免了仅使用来自单个 ICESat-2 或 GEDI 轨道的高程值的情况。为了达到更真实的视觉景观表示，我们仅将源自原始 CopernicusDEM 的表面粗糙度添加到插值后的地形值中。粗糙度或地形位置指数是像素高程与其八个相邻像素平均高程之间的差值。由于这些值有时可能达到数米，我们将其限制在低分辨率 DTM 的初始阈值（用于 PMF 滤波器）的范围内。因此，如果初始阈值为 1.2 米，则范围为 -0.6 至 0.6 米，所有大于 0.6 米的 TPI 值将被设置为 0.6 米，所有小于 -0.6 米的值将被设置为 -0.6 米。在最坏的情况下，这会向 DEM 添加随机噪声，类似于未插值的 CopernicusDEM 高程值中存在的噪声。然而，在最佳情况下，它能够反映真实的地形特征，例如树冠下的沟渠或小运河。总体而言，这些新增区域较小且分布均衡（大致呈零均值），不会影响精度。

精度验证

使用澳大利亚、佛罗里达州、印度尼西亚、拉脱维亚、马绍尔群岛、墨西哥、荷兰、波兰和英国的公开本地机载激光雷达参考数据集对 DeltaDTM 数据集进行验证。这些数据集总面积达 78106 平方公里，^{覆盖了全球平均海平面附近或以下的沿海地区。所有数据集均使用 GDAL}³⁶从其本地垂直参考点垂直重投影到平均海平面（使用 EGM2008 大地水准面）。

DeltaDTM 在所有土地覆盖类型综合应用中表现最佳，偏差为 0.01 米，平均绝对误差 (MAE) 为 0.45 米，均方根误差 (RMSE) 为 0.74 米。DeltaDTM 91% 的数据点与参考表面的偏差在 1 米以内，98% 在 2 米以内，100% 在 5 米以内。DiluviumDEM 的表现次之，FABDEM 紧随其后——CoastalDEM 与之非常接近，但在 1 米以内的偏差比例上略逊一筹。DiluviumDEM 的偏差为 0.02 米，MAE 为 0.72 米，RMSE 为 1.22 米，79% 的数据点与参考表面的偏差在 1 米以内。FABDEM 的偏差为 0.65 米，MAE 为 1.05 米，RMSE 为 1.96 米，71% 的数据点与参考表面的偏差在 1 米以内。

引用：

Pronk, M., Hooijer, A., Eilander, D. et al. DeltaDTM: A global coastal digital terrain model. Sci Data 11, 273 (2024).https://doi.org/10.1038/s41597-024-03091-9

数据集引用：

Pronk, Maarten (2024): DeltaDTM: A global coastal digital terrain model. Version 3. 4TU.ResearchData. dataset.https://doi.org/10.4121/21997565.v3

—

GEE部分下载代码

var delta_dtm = ee.Image("projects/sat-io/open-datasets/DELTARES/deltadtm_v1"),    roi = /* color: #98ff00 */ee.Geometry.Polygon(        [[[115.52822350926307, 40.0223076464916],          [117.02236413426307, 28.356303383109573],          [128.79970788426306, 31.684884862973604],          [127.12978600926307, 37.973382302325945],          [123.52627038426307, 40.491822280168726]]]);var elevation = delta_dtm.select('b1');elevation = elevation.updateMask(elevation.neq(10));//设置底图var snazzy = require("users/aazuspan/snazzy:styles");snazzy.addStyle("https://snazzymaps.com/style/132/light-gray", "Grayscale");var elevationVis = {  min: 0,  max: 10.0,  // cmocean deep  palette: ["281a2c", "3f396c", "3e6495", "488e9e", "5dbaa4", "a5dfa7", "fdfecc"]};Map.setCenter(117.07, 32.76, 7);  // Map.addLayer(elevation, elevationVis, 'DeltaDTM');Export.image.toDrive({  image: elevation.clip(roi),//选择导出影像的波段  description: 'dem',//选择导出云盘的文件夹名称  crs: "EPSG:4326",//坐标系  scale: 30,//空间分辨率  region: roi,//研究区  maxPixels: 1e13,//最大像元个数  folder: 'dem'});