数字高程模型（Digital Elevation Model，DEM）作为最基础的地理信息产品之一^[1]，在基础工程建设、灾害监测、资源调查等领域发挥着重要作用^[2]。传统光学遥感和干涉合成孔径雷达（Interferometric Synthetic Aperture Radar，InSAR）是目前全球尺度DEM数据获取的重要手段^[3]，然而由于对地观测技术的局限性，获取的全球尺度DEM数据中不可避免地包含植被、建筑物等信息，这将影响DEM对实际地形的表达精度。科学研究中所使用数据的精度将直接影响研究结果的可靠性，因此充分了解所使用的全球尺度DEM数据的误差情况并利用可靠手段提高其精度至关重要。

ASTER GDEM是一款通过光学立体摄影测量得到的全球尺度DEM产品，自ASTER GDEM发布以来，许多学者都聚焦于对其精度的评价。目前常用的精度评价方法是选择高精度的DTM或GNSS测量点作为参考数据，然而这两种数据采集成本高、工作量大。ICESat（Ice, Cloud, and land Elevation Satellite）卫星自2003年发射以来获取了大量全球尺度的高精度高程数据，ICESat-2于2018年发射，它是对ICESat工作的延续与改进， ICESat-2获取的数据在数据密度和精度相较于ICESat都得到较大提升^[4]。ICESat-2发射同年，美国发射了全球生态系统动力学调查雷达（GEDI），其主要用于全球尺度的森林冠层高度、垂直结构以及地表高程的精准测量^[5-6]。二者能够在大范围区域内迅速地获取大量的高精度高程点^[7]。研究表明，ICESat-2数据产品的地表高程精度优于0.8 m^[8-9]，GEDI数据产品的地表高程精度优于1.8 m^[10-11]，精度远远高于ASTER GDEM。这些星载激光雷达高程数据为DEM精度评价提供了一种新途径。

随着DEM应用的日益广泛，DEM误差修正成为学者们关注的研究热点。如Gallant等^[12]利用Landsat TM数据识别出植被覆盖地块，通过估算每个地块的植被高从而移除SRTM（Shuttle Radar Topography Mission） DEM中的植被高信号；Baugh等^[13]利用ICESat数据在林区的测量结果估算SAR（Synthetic Aperture Radar）信号对植被的穿透深度从而移除SRTM DEM中的植被高信号；Tan等^[14]以GEDI冠层高度为参考，利用InSAR相干散射模式估算森林高度，在此基础上建立了基于森林高度和植被覆盖度的线性回归模型，并利用该模型移除TanDEM-X（TerraSAR-X add-on for Digital Elevation Measurements）DEM中的植被高信号；张晨等^[15]基于ICESat-2数据使用随机森林算法反演得到植被高从而扣除TanDEM-X DEM中的植被高信号。上述研究主要是针对DEM误差中的植被高信号，没有考虑其他地物信号及地形地貌因素，这限制了方法在多种土地覆盖类型下的应用。考虑各种影响因素的线性回归方法也常用于DEM误差修正。如杜小平等^[16]通过建立ICESat数据与ASTER GDEM的一元线性回归模型实现对DEM的修正；Su等^[17]通过建立DEM误差随树高、植被覆盖率和地形坡度变化的线性回归模型修正SRTM DEM。此外，秦臣臣等^[18]采用反距离加权、普通克里金、地形转栅格和自然邻域4种空间插值方法构建误差曲面从而修正SRTM DEM和ASTER GDEM。然而，单纯地使用空间插值或线性回归模型并不足以挖掘各种影响因素与DEM误差的关系从而限制了误差修正的效果。

XGBoost机器学习模型于2016年由Chen等^[19]提出，该模型是在梯度提升树模型（Gradient Boosting Decision Tree）的基础上进行优化而得到的一种高效模型，强大的数据拟合能力使得它能够把大量非线性关系的特征要素拟合在一起，具有精度高、模型简单、支持并行化以及能有效防止过拟合等优势，并且借鉴了随机森林算法，支持列抽样，训练速度快。单个星载激光雷达足迹通常沿地面轨道均匀分布而不是随机分布，这会产生强烈的条带效应，结合ICESat-2和GEDI两种不同的星载激光雷达数据不仅可以缓解条带效应还可以提高高精度参考高程点的地面覆盖程度及数量^[6]。为此，本文以ICESat-2及GEDI数据为参考，评价ASTER GDEM的精度，并分析了土地覆盖类型、高程、坡度、起伏度及植被覆盖率与ASTER GDEM精度的关系。在此基础上，提出了一种结合XGBoost与空间插值的DEM修正方法。

本文以福建泰宁为研究区对ASTER GDEM进行精度评价和修正。福建泰宁位于我国东南部地区武夷山脉中段的东南侧，26.5°N—27.5°N，116.5°E—117.5°E，占地总面积1539.38 km²（图1）。

泰宁涵盖多种土地覆盖类型，地形以山地和丘陵为主，海拔差异较大，山地多平原少，河流众多，水系发达，植被以中南亚热带常绿阔叶林为主，山地占38.7%，丘陵占59.9%，河谷盆地占1.4%，复杂的地形地貌和茂密的植被覆盖为本文的研究提供了一个良好的实验案例。

ASTER GDEM是由日本和美国联合研制并免费公开的全球数字高程产品，该产品是通过光学立体摄影测量的手段获取，由于包含地表建筑物及植被信息，严格意义上ASTER GDEM属于DSM产品^[20]。ASTER GDEM有v1、v2、v3共3个版本，覆盖了全球83°S—83°N范围内的所有陆地区域，分辨率为1″（约30 m），数据的获取时间均为2000—2009年。GDEM v1是在无控制点的情况下利用全球范围光学立体像对全自动生产完成，GDEM v2较v1采用了更加先进的算法，提高了高程精度和空间分辨率精度，GDEM v3在v2的基础上新增了大量光学立体像对，覆盖范围及高程精度得到一定提升。本文采用ASTER GDEM v3数据^[21]（https://earthdata.nasa.gov/），垂直标称精度为±20 m，参考椭球体为WGS84。

ICESat-2卫星搭载了光子计数激光测高仪（Advanced Topographic Laser Altimeter System，ATLAS），ATLAS首次在卫星平台应用了微脉冲多波束光子计数激光雷达技术，新技术的应用使得ICESat-2在林下地形及冠层高度测绘方面具有巨大优势^[4,22]。每次脉冲ATLAS产生6个独立光束，每条独立光束可以沿轨道获取间隔约0.7 m、直径约 17 m 的重叠光斑。6个独立光束沿轨排列成平行的3对进行对地观测，每对又分为能量功率比为4：1的强弱光束。ICESat-2/ATLAS包括22种标准数据产品，本文主要采用ATL03及ATL08数据产品^[23-24]，ATL03 提供了接收光子的经纬度、椭球高度（以WGS84为参考椭球体）及光子分类等基础信息，是生成其他产品的基础数据，而ATL08数据是在ATL03数据的基础上处理得到的每100m统计单元的陆地及植被参数（https://nsidc.org/data/icesat-2/）。

GEDI搭载了一种全波形多光束激光高度计（Full-waveform, Multibeam Laser Altimeter），该高度计共携带3个激光发射器，产生沿轨共8条平行光束（4条覆盖光束、4条全功率光束）^[10]。每条独立光束可以沿轨道获取间隔约60 m、直径约25 m 的光斑。GEDI的数据产品分为4个级别，本文主要采用GEDI的二级数据产品L2A^[25],该产品提供了沿轨60 m每个光斑中心位置的地表高程和冠层分位数高度指数^[26]（https://lpdaac.usgs.gov/products/gedi02_av001/）。

本文获取了覆盖研究区全域的ICESat-2及GEDI数据（图1），共包含研究区12个条带的ATL03及ATL08数据产品及9个条带的L2A数据产品。

30 m全球地表覆盖数据（GlobeLand30）是我国自主研发的30 m空间分辨率全球土地覆盖数据^[27]。Globleland30共有3个版本，对应不同年份的全球土地覆盖类型。GlobeLand30数据研制所使用的分类影像主要是30 m多光谱影像，包括Landsat-5专题制图仪、Landsat-7增强型专题制图仪、OLI多光谱影像和中国环境减灾卫星（HJ-1）多光谱影像，2020版数据还使用了16 m分辨率高分一号（GF-1）多光谱影像。GlobeLand30数据共包括10个一级类型，分别是：耕地、林地、草地、灌木地、湿地、水域、苔原、人造地表、裸地、冰川和永久积雪，总体分类精度为85.72%，本文使用的是Globleland30的2020版本（图2(a)）。

植被连续场（Vegetation Continuous Field，VCF）提供地表固定空间分辨率内的植被覆盖率数据^[28]。为匹配ASTER GDEM的空间分辨率和时间尺度， 本文使用的VCF是由2000—2009年的Landsat VCF经过中值合成得到（图2(b)）。Landsat VCF由Landsat-5专题制图仪和Landsat-7增强型专题制图仪获取的所有波段数据处理得到，包含30 m×30 m像素内高度大于5 m植被的植被覆盖率^[29-30]。

坡度和起伏度使用ArcGIS软件对ASTER GDEM进行处理得到，其中起伏度是通过计算ASTER GDEM九宫格像素内最大高程值与最小高程值的差值获得。栅格数据（ASTER GDEM、GlobeLand30、Landsat VCF、坡度、起伏度）的空间分辨率是30 m，存在着空间坐标系不同以及像素位置不匹配的可能性。为此，本文将上述栅格数据转换到ASTER GDEM空间坐标系中，并以ASTER GDEM像素为准采用双线性插值方法匹配所有栅格数据的像素位置。

ATL08数据产品提供了ICESat-2卫星沿轨100 m采样间隔的地表高程值，而ASTER GDEM的空间分辨率为30 m。为匹配二者的分辨率，使用开源工具Phoreal（Photon Research and Engineering Analysis Library， https://github.com/icesat-2UT/PhoREAL）将ATL08光子分类索引与ATL03数据关联并对ATL03光子进行标记从而建立ATL08与ATL03间的联系，通过计算ATL03数据沿轨每30 m步长范围内地表光子的高程平均值对ATL03光子进行重采样，得到ICESat-2沿轨每30 m步长中心点的地表高程，这种平均值更能反应沿轨每30 m内的地 表高程以匹配ASTER GDEM的空间分辨率。 GEDI_L2A数据集中提供的quality_flag和 degrade_flag 分别表示数据质量标志和卫星指向及定位状态降级标志，通过筛选出quality_flag参数等于1以及 degrade_flag参数等于0的光斑得到高质量GEDI高程点。最终得到了研究区内ICESat-2及GEDI高程点共计34 269个。为了消除噪声高程点的影响，首先根据Huang等^[31]的方法删除高差绝对值大于100 m的高程点，再根据拉依达准则进一步删除高程点。此外，由于2种星载激光雷达高程点在研究区有交叠，对落到同一个像素中的多个高程点取平均值。最终剩余33 882个高程样本点。本文计算了每个ICESat-2及GEDI高程点与ASTER GDEM高程的差值。

通过统计分析所有高程点位置处的DEM误差评价原始ASTER GDEM的整体高程精度。此外，本文提取了所有高程点位置的土地覆盖类型、高程、坡度、起伏度及植被覆盖率，通过分带误差统计分析ASTER GDEM精度与上述影响因素的关系。

为了有效地评价修正方法的效果，本文随机 将ICESat-2及GEDI高程点分成75%、20%及5% 3份数据集，其中75%的高程点用于建立研究区的DEM绝对误差曲面及DEM误差正负权重曲面，5%的高程点用于建立多元非线性回归方程，20%的高程点将作为检验数据集用于评价修正后ASTER GDEM的高程精度。

本文采用平均误差（ME）、平均绝对误差（MAE）、均方差（RMSE）以及决定系数（R²） 4种分析指标评价原始ASTER GDEM的精度及DEM修正方法的效果，其中前3个指标是广泛使用的参数，其定义在此不再赘述。决定系数（R²）也称为拟合优度，在文中仅用于修正后的数据精度评价，它反映了2种数据的接近程度^[32]，R²越接近1，2种数据的相关性越高，计算公式如下：

式中：Hi,ICESat-2&GEDI表示第i个高程点的地表高程值；Hi,Corrected ASTER GDEM表示修正后的ASTER GDEM在第i个高程点对应位置的高程值；H-i,ICESat-2&GEDI表示检验数据集中所有高程点的平均值；n表示检验数据集中高程点的总数。

地表覆盖物及地形地貌是研究区ASTER GDEM误差的主要来源。本文基于高程、坡度、起伏度、植被覆盖率及地理位置（经纬度）这5种因素与DEM误差之间的关系建模，来实现对DEM的修正。在上述各种影响因素的交互影响下，DEM在不同环境下的误差正负较为随机，特别在样本点较少的小范围研究区，很难寻找上述影响因素与误差正负的相关性。本文利用Spearman相关系数评估了上述影响因素与DEM误差及DEM绝对误差的相关性，Spearman相关系数是衡量2个变量的依赖性的非参数指标，它利用单调方程评价2个统计变量的相关性，计算公式如下：

式中：xi为第i个高程样本点位置处的精度影响因素（如高程）数值；yi为第i个高程样本点位置处的DEM误差或DEM绝对误差；x-表示所有高程样本点位置处精度影响因素数值的平均值；y-表示所有高程样本点位置处DEM误差或DEM绝对误差的平均值。

从计算的Spearman相关系数（表1）可见，除了经纬度之外，其他因素与DEM绝对误差的相关性均高于其与DEM误差的相关性，因此本文选择将DEM绝对误差作为XGBoost模型的目标值，高程、坡度、起伏度、植被覆盖率及地理位置（经纬度）作为输入值，借助XGBoost模型对研究区的DEM绝对误差进行回归预测，从而得到研究区的DEM绝对误差曲面。

选择DEM绝对误差作为XGBoost模型的目标值，其中一个关键性问题即如何确定误差的正负。根据地理学第一定律，误差为负的高程样本点附近DEM误差大概率为负值，同理误差为正的高程样本点附近的DEM误差大概率为正值。基于上述思想，本文利用一定数量的高程样本点插值出研究区的DEM误差正负权重曲面，首先将正误差样本点的DEM误差值标准化到范围0~1，负误差样本点的DEM误差值标准化到-1~0，再根据每个样本点经过标准化后的DEM误差正负权重利用反距离加权法插值出研究区的DEM误差正负权重曲面，这个权重曲面将用于有效评估研究区DEM误差正负情况。

通过建立DEM绝对误差预测值及DEM正负误差权重值与DEM误差值之间的多元非线性回归方程得到研究区DEM误差曲面，叠加原始ASTER GDEM和DEM误差曲面消除DEM误差实现DEM修正。图3为ASTER GDEM修正及精度评价流程。

本文以ICESat-2及GEDI高程点作为参考数据对ASTER GDEM高程误差进行统计分析（图4），误差通过计算ICESat-2及GEDI高程点与ASTER GDEM高程值的差值求得。ASTER GDEM的误差整体呈正态分布，平均误差为-3.463 m，平均绝对误差为12.930 m，均方差为16.695 m。

不同土地覆盖类型下的ASTER GDEM误差统计情况表明（表2）：不同土地覆盖类型下的ASTER GDEM高程精度存在显著差异，不同土地覆盖类型下的MAE及RMSE由大到小依次是：林地、草地、耕地、水域、人造地表和湿地，这种精度差异产生的原因主要是不同土地覆盖类型的地表覆盖物及地形地貌。除人造地表外，其他土地覆盖类型下的DEM平均误差均为负值，这是由于ASTER GDEM光学摄影测量的数据采集手段对地表地物没有穿透性导致ASTER GDEM包含了地表各种地物信息。人造地表的ME为正的原因是ASTER GDEM与星载激光雷达数据的采集时间不同，这段时间内人造地表区域内建筑物的变化导致的。

基于高程、坡度、起伏度及土地覆盖率（VCF）对ASTER GDEM进行分带精度统计情况如图5所示。从坡度分带及起伏度分带对应的DEM平均误差及平均绝对误差的变化情况可以看出，ME随着坡度及起伏度的增加而减小，MAE随着坡度及起伏度的增加而增大，这说明随着坡度及起伏度的增加，DEM误差越来越大且越来越高估地表高程。

根据VCF分带的DEM平均误差及平均绝对误差的变化情况可看出，MAE随着VCF的增加而增大，ME随着VCF的增加先增大后减小再增大。ME从VCF 0~10%至10%~20%迅速增大，之后随着VCF的增加缓慢减小，但当VCF大于60%~70%后又开始增大。

根据高程分带的DEM平均误差及平均绝对误差的变化情况可以看出，MAE从高程200~300 m至300~400 m迅速增大，从高程300~400 至600~700 m缓慢增大，从高程600~700 m至800~900 m基本不变，但从高程大于800~900 m后又开始增大。ME随着高程的增加先减小后增大，增大的临界值在高程700~800 m之间，这是由于随着高程的增加地表地形地貌及植被覆盖复杂度增加导致DEM越来越高估地表高程，但当高程达到一个阈值之后，温度气压的降低导致地表植被覆盖率的下降，这使得植被对DEM高估地表高程的影响降低。

利用XGBoost模型回归预测得到的DEM绝对误差曲面在空间分布上有明显的条带性（图6）。通过对检验数据集高程点提取得到的DEM绝对误差曲面值与实际DEM绝对误差的差值进行统计分析，评价XGBoost对绝对误差的回归预测精度（表3）。

利用多元非线性回归的方式建立本文3.2节提出的方法插值得到的DEM误差正负权重及DEM绝对误差与DEM实际误差之间的方程：

式中：x为每个对应位置处的DEM绝对误差曲面值；y为每个对应位置处的DEM误差正负权重曲面值；σz为每个对应位置处最终回归得到的DEM修正值。

利用上述方程得到DEM误差曲面，通过叠加DEM误差曲面与原始ASTER GDEM得到修正后的DEM。为了更直观地展示修正前后的DEM 变化，本文对修正前后的DEM进行了剖面分析 （图7）。图7是ICESat-2经过研究区的一条轨迹的剖面，叠加了修正前后的ASTER GDEM。其中不同颜色的点表示不同分类标签的ICESat-2光子点（地面（棕色）、冠层（绿色）和噪声（灰色），紫色和红色数据分别为ASTER GDEM未修正高程和修正高程），可以看出修正后的ASTER GDEM更加贴近真实地表，对地表高程的低估和高估都得到一定修正。在图7所示样带轨迹中，原始ASTER GDEM经过修正后的平均误差由-1.152 m降低到了-0.047 m，平均绝对误差由8.707 m降低到了3.994 m，均方差由10.769 m降低到了5.383 m。

图8（a）为ASTER GDEM在修正前后的整体误差分布情况。修正前，ASTER GDEM的平均误差为-3.463 m，说明ASTER GEDM在研究区有较大的负偏差且普遍高估地表高程。修正后，ASTER GDEM的平均误差降低到-0.233 m，负偏差显著降低，平均绝对误差和均方差分别由12.930 m和16.695 m降低到了9.564 m和12.762 m，数据精度得到了有效改善。此外，修正后DEM误差更加集中地分布在0 m附近，说明地表覆盖物对DEM的影响得到降低。图8（b）和图8（c）为ASTER GDEM在修正前后的精度验证散点图，R²由0.995提高到了0.999。

修正前后不同土地覆盖类型下的ASTER GDEM精度变化情况如表4所示。修正前，ASTER GDEM的平均误差在林地最大（-4.491 m），在人 造地表最小（0.522 m）。经过修正后，耕地、林地、 草地、湿地、水域及人造地表的平均误差都得到 显著降低，分别由-0.435、-4.344、-2.464、-4.232、 -3.560、0.738 m降低到了-0.333、-0.349、0.956、 -1.147、-0.223及-0.124 m。修正前，林地ASTER GDEM的平均绝对误差和均方差最大，分别为14.077 m和17.866 m，湿地ASTER GDEM的平均绝对误差和均方差最小，分别为6.275 m和7.728 m。经过修正后，不同土地覆盖类型下的ASTER GDEM平均绝对误差及均方差都得到一定改善，耕地、林地、草地、湿地、水域及人造地表的平均绝对误差分别降低了29.96%、24.73%、20.79%、58.02%、62.41%、31.47%，均方差分别降低了27.30%、23.25%、17.69%、43.96%、53.02%、30.41%。消除地表覆盖物及地形地貌对DEM的影响是DEM修正的本质，因此耕地、林地、草地、人造地表这种地表覆盖物复杂且地形地貌变化大的地方修正效果较差，湿地和水域这种地表覆盖物简单且地形平坦的地方修正效果较好。

（1）已有的DEM修正方法多采用线性回归模型或空间插值并且多数研究仅考虑了植被高信号，没有将地形因子及土地覆盖类型综合考虑进去。本文提出的DEM修正方法不仅考虑了多种影响因素，并且结合了XGBoost和空间插值，能够更深入地挖掘DEM误差与多种影响因素的关系。此外，已有研究使用的数据多为单一类型的星载激光雷达且以ICESat为主，而本研究结合了2种不同类型的星载激光雷达数据，有效提高了地表高程参考点的空间覆盖范围和数量。

（2）本研究目前涉及的实验区域面积相对较小，地形地貌和植被覆盖特征样本数量可能存在不足，本文提出方法的普适性和可推广性尚不明确（如：修正方法能否直接应用或移植到其他类似的区域）。为此，后续我们将在范围更大的区域中开展研究，并利用实地采集的高精度地形数据（如：通过地面和无人机激光雷达获取的高密度点云）对修正方法的效果做进一步的验证分析。

（3）我国于2022年8月4日成功发射首颗陆地生态系统碳监测卫星“句芒号”，搭载包括了星载激光雷达在内的多种遥感传感器，本文提出的DEM精度评价及修正方法，将来也可以应用在国产星载激光雷达生产的数据产品上。

本研究基于ICESat-2及GEDI评估了ASTER GDEM v3的高程精度并分析了其精度与土地覆盖类型、高程、坡度、起伏度及植被覆盖度VCF之间的关系。在此基础上提出一种结合XGBoost和空间插值的ASTER GDEM修正方法。结果分析表明：

（1）本文提出的DEM精度评价及修正方法使用了新型星载激光雷达ICESat-2及GEDI数据，对地形误差与多种影响因素之间的非线性关系进行拟合、建模，得到了较为精准的修正结果。

（2）研究区原始ASTER GDEM高程误差整体呈正态分布，平均误差为-3.463 m，平均绝对误 差和均方差分别为12.930 m和16.695 m。原始 ASTER GDEM高程精度随着高程、坡度、起伏度及VCF的增大呈现下降趋势，且不同土地覆盖类型下的高程精度存在显著差异。

（3）本文提出的DEM修正方法实现了对ASTER GDEM精度的有效改善。修正后，ASTER GDEM的平均误差由-3.463m降低到了-0.233 m。DEM的整体平均绝对误差降低到了9.564 m，修正了26.04%，整体均方差降低到了12.762 m，修正了23.56%。耕地、林地、草地、湿地、水域及人造地表的DEM平均绝对误差分别降低了29.96%、24.73%、20.79%、58.02%、62.41%、31.47%，DEM均方差分别降低了27.30%、23.25%、17.69%、43.96%、53.02%、30.41%。