数字高程模型(digital elevation model，DEM)是地表起伏形态的数字化表达，也是进行各类地学分析的重要基础数据，已被广泛用于水文、地质、气象及军事等各个领域^[1-6]。伴随着对地观测技术的不断发展，获取DEM数据的能力日益增强，各研究领域对高质量DEM的需求也随之增长^{[4, 7-11]}。目前，由于观测手段、处理方式和地形条件等因素的影响，DEM仍然存在异常值、噪声、空洞等质量问题，直接影响数据的应用潜力^[9]。多源DEM融合方法通过综合不同数据之间的互补信息，能获取更准确、全面、可靠的DEM，实现现有数据集的质量提升。

规则格网DEM因数据形式简单、处理分析便捷而成为最常用的数据形式。在全球尺度的栅格DEM产品中，ASTER GDEM(advanced spaceborne thermal emission and reflection radiometer global digital elevation model)和SRTM DEM(shuttle radar topography mission digital elevation model)的应用最为普遍^{[8, 12]}。ASTER GDEM是根据Terra卫星先进星载热发射和反辐射计的立体像对而生成的，以高分辨率(30 m)覆盖全球约99%的陆地表面。在发布的2个版本中，ASTER GDEM v2(第2版数据)的分辨率和精度相对较高。但由于光学成像易受云雾等天气的影响，数据中仍然存在大量异常值和噪声(图 1(b))，使其精度饱受争议^[13]。SRTM DEM是根据美国奋进号航天飞机获取的雷达影像所制作的高程产品，覆盖全球80%以上的陆地表面，因其精度在平坦地区稳定而成为最常用的DEM数据源^[14-15]。然而受限于雷达侧视成像模式，SRTM DEM易在坡度较大的区域形成数据空洞。USGS根据不同处理级别发布了90 m分辨率的SRTM3和30 m分辨率的SRTM1。目前，SRTM3 v4.1是使用最广泛的DEM，在SRTM v2.1的基础上结合Reuter等的算法填充空洞所得，然而这版数据的空洞填充质量并不稳定，在山区的垂直精度差异较大^{[13, 16]}。SRTM1从2014年9月起面向全球发布，具有较高和较稳定的垂直精度，但仍有空隙和异常(图 1(a))。

除了栅格DEM数据集之外，全球尺度还有丰富的矢量高程产品。例如，美国国家航空航天管理局(NASA)在2003年发射了冰、云和陆地高程卫星ICESat (ice，cloud，and land elevation satellite)，搭载的地学激光测高系统GLAS(geoscience laser altimeter system)在2003-2009年获取了全球范围的高程点(图 1(c))，其垂直精度在平坦地区可达厘米级^[17]。

目前，针对DEM产品的质量问题，国内外学者开展了大量融合处理工作。一方面，对于噪声和异常值，主要有以下研究：文献[18]基于插值算法，利用ICESat GLAS激光测高数据得到差值面，尝试校正ASTER GDEM数据，然而，由于ICESat GLAS数据分布稀疏，在离采样点较远的位置误差较大；文献[10]通过在频率域融合低分辨率DEM的高频成分和高分辨率DEM的低频成分实现去噪，但结果存在细节信息丢失现象。另一方面，针对DEM的空洞问题，现有研究工作主要包括：文献[12]提出TIN差分曲面，解决了空洞填充过程中不同数据间的高程偏差；文献[19]针对SRTM及辅助数据在空洞周边区域的特点选择相应的插值方法，提出多源DEM融合策略，得到了无缝DEM产品SRTM v4.1；文献[20-21]构建高精度曲面建模方法，支持融入多源数据，可处理DEM的空洞问题。此外，也有学者同时针对空洞和噪声发展融合方法：文献[8]利用数据拼接和平滑算法，通过融合多个DEM数据集生成“EarthEnv-DEM90”产品，但未充分考虑输入数据集的精度差异；文献[22]提出一种基于正则化变分的多尺度融合方法，能够同时克服噪声和空洞的影响，有效提升数据分辨率，但融合结果的分辨率和精度水平受限于输入DEM。

总的来说，现有的融合方法没有充分考虑多源数据之间在获取方式、分辨率、覆盖范围及精度上的差异，难以同时处理噪声、异常值、空洞、分辨率差异等问题^{[8, 10, 12, 19]}。此外，目前的融合方法也较少引入高精度矢量数据，在一定程度上限制了融合产品的质量。对于多源DEM融合而言，既要考虑形式和精度上的差异，又要综合各数据集的优势，才能得到精度更高的DEM产品。因此，结合SRTM1、ASTER GDEM v2及ICESat GLAS数据的特性，本文将在笔者前期方法的基础上^[9]，进一步提出基于高差拟合神经网络的多源DEM融合方法，改善模型的自适应特性，得到空间无缝DEM。

1 数据融合方法 1.1 多源DEM数据融合框架

本文以SRTM1为基础数据，利用ASTER GDEM v2填补其数据空洞。虽然ASTER GDEM v2覆盖范围较广，但有异常值和噪声。为了充分发挥ASTER GDEM v2在覆盖范围上的优势，利用ICESat GLAS校正其垂直精度。ICESat GLAS的数据沿轨道点状分布，间距大，传统插值方法误差较大。因此，本文以点面融合为基本思想，利用神经网络训练点状ICESat GLAS与面状ASTER GDEM v2的误差统计关系，实现对后者的精度校正。在校正过程中，根据DEM的空洞率(空洞高程点占所有高程点的比例)采取不同的处理策略。对于空洞率小于0.05%的SRTM1数据，缺失范围很小，直接利用基于TIN差分曲面的算法填充空洞。而对于空洞率大于0.05%的数据，则先经过精度校正，再进行空洞填补。根据试验分析，ICESat GLAS与ASTER GDEM v2的值都表示为高程，可能存在一定的线性相关性。因此，为了使网络能更好地学习数据关系，利用ICESat GLAS与ASTER GDEM v2高差进行训练，以削弱数据间的相关性。

综上所述，本文方法的整体流程分为以下3个步骤(如图 2所示)：①数据预处理。实现多源DEM的坐标统一，并滤除ICESat GLAS的异常值。②高差拟合神经网络点面校正。在SRTM1的空洞范围内，利用ICESat GLAS对ASTER GDEM v2进行精度校正。③TIN差分曲面无缝填充。在校正后的ASTER GDEM v2基础上，以TIN差分曲面方法填补SRTM1空洞，并对残留的异常值进行滤波，得到最终融合结果。

1.2 数据预处理

预处理阶段主要是ICESat GLAS数据坐标系转换和异常值滤除。具体为：①进行坐标转换。SRTM1和ASTER GDEM v2采用的是WGS84/EGM96坐标系，而ICESat GLAS采用TOPEX/Poseidon坐标系。为了消除不同高程系统带来的系统偏差，需要统一坐标系。由于两种坐标系在经纬度方向上的水平偏差不到1 m，只在垂直方向上进行ICESat GLAS的高程转换^[23]。利用式(1)将椭球高h_TOPEX转换成椭球高h_WGS-84，再转换为大地水准高H。②进行异常值滤除。结合常用的参数标准筛选出ICESat GLAS后^[18]，再利用式(2)的坡度自适应思想设置高程阈值，滤除受遮挡物影响的异常值^[24]，确保对于SRTM1或ASTER GDEM v2的高程绝对误差在合理范围内。

式中，h_TOPEX、h_WGS-84分别为对应于TOPEX、WGS-84的椭球高；0.707为两椭球体的垂直偏差；H为对应于EGM96模型的大地水准高；N为大地水准面差距；单位均为米。

式中，S为坡度向量；i为相应高程点位置；y′_v为参考高程，即SRTM1或ASTER GDEM v2高程；y_u为ICESat GLAS高程；T为设置的阈值。

1.3 基于高差拟合神经网络的点面校正方法

人工神经网络(artificial neural network，ANN)只需确定非常量性参数，便能自组织、自学习、自适应地挖掘数据统计特征。因其强大的拟合能力、学习能力和较强的鲁棒性，被广泛地应用于回归问题中^[25-26]。本文通过ASTER GDEM v2和ICESat GLAS构建训练样本集，挖掘二者的统计关系，得到ASTER GDEM v2的高程预测值，实现其精度校正，具体流程如图 3所示。将ASTER GDEM v2的经纬度坐标和地形坡度信息作为网络输入、ICESat GLAS数据点与对应ASTER GDEM v2的高程差值作为网络输出，将网络输出与ASTER GDEM v2相加，即为校正后的高程值。

顾及ICESat GLAS和ASTER GDEM v2的水平位移差及异常值的影响，通过邻域窗口选取样本来构建输入矩阵，经过归一化处理后，将输入向量和目标向量输入神经网络。整个模型由n个输入层神经元、s个隐含层神经元和p个输出层神经元构成，通过权值参数和激活函数连接输入向量和目标向量，实现信息传递和学习。模型中，设输入层神经元为x=(x₁, x₂, …, x_n)，用于输入数据的接收；隐含层和输出层神经元分别为h=(h₁, h₂, …, h_s)，o=(o₁, o₂, …, o_p)，用于接收数据的加工和处理(式(3))。由于输入层和输出层之间存在一定的映射关系，基于梯度下降的误差逆传播算法，将输出层误差通过逐层回传来调整神经元的连接权和阈值^[27]，使输出值o与真实值y的误差E达到最小(式(4))。

式中，w_ij^h、w_jk^o分别为输入层神经元x_i与隐含层神经元h_j之间的连接权、隐含层神经元h_j与输出层神经元o_k的连接权；β_j、τ_k为对应的阈值；f₁(·)、f₂(·)分别为隐含层、输出层的激活函数，反映神经元输入和输出之间的函数关系。

式中，E为模型误差；n为迭代次数；η为学习速率。其中，学习率关系着连接权的更新，其值越大，往往收敛速率越快，但易造成震荡效应。

模型参数设置如下：邻域窗口为3×3；隐含层神经元个数为6；隐含层、输出层的激活函数分别为“tansig”和“purelin”函数；学习率为0.001；输入层和输出层的神经元个数根据变量数而定。

1.4 基于TIN差分曲面的无缝填充方法

在数据融合过程中，多源DEM间的垂直偏差会造成边缘出现“断裂”现象，因此，利用文献[12]提出的基于TIN差分曲面的算法来处理这种偏差。此算法的基本思想是，在空洞区域内，假设校正后的ASTER GDEM v2与SRTM1的地形起伏保持一致，由ASTER GDEM v2获取的差分曲面同样可以用于SRTM1。因此，本文在SRTM1空洞区域周边提取宽度为5个高程点的缓冲区，构建空洞区域基面，进而通过ASTER GDEM v2基面和校正后的ASTER GDEM v2得到差分曲面，应用于SRTM1，过程表示如下

式中，S_o、A_o分别为SRTM1、ASTER GDEM v2数据的基面；A_C为空洞范围内校正后的ASTER GDEM v2的高程；S_C为最终填充SRTM1的曲面。鉴于填充数据可能存在噪声和异常值，在填充前对曲面进行自适应平滑滤波，且在填充后对连接区域进行轻微羽化，实现接缝处的平滑过渡。

2 试验结果与分析 2.1 模型验证

对神经网络校正模型的空间预测能力进行验证。本文选用了交叉验证^[9]和中心块验证两种方式测试模型精度。交叉验证将ASTER GDEM v2覆盖的有效ICESat GLAS分成10等份，以其中的9份作为训练数据，剩下的1份则作为测试数据，循环试验直到所有的数据都曾作为测试数据。由于ICESat GLAS数据在沿轨道方向和轨道间分布差异较大，则同时采用中心块验证的方法(如图 4所示)，将训练和测试数据分为空间上互不重合的部分，进一步验证模型的可靠性和稳定性。具体地，将ASTER GDEM v2中心区域的ICESat GLAS作为测试数据，其余的ICESat GLAS作为训练数据。表 1-2分别给出了交叉验证和中心块验证的精度评价情况，包括原始ASTER GDEM v2、多元线性回归方法(multiple linear regression，MLR)及本文方法校正ASTER GDEM v2的结果，分别表示为ASTER、MLR、ANN。由于均方根误差(root mean square error，RMSE)对异常值敏感，能很好地反映数据精度，将其作为定量评价指标(式(6))。

式中，s_a、z_a分别为第a个参考数据ICESat GLAS及待评价数据集对应点的高程值；N为测试高程点数。

从表 1、表 2均可看出，在所有验证试验中，神经网络模型结果相对于原始数据和MLR结果更优，表明其能较好地构建数据之间的误差关系。表 2的中心块验证，因训练数据与测试数据在空间覆盖范围差异较大，能更好地验证模型的效果，且神经网络模型结果优势更为突出，有利地验证了其空间预测能力更强。

为了进一步说明精度校正的必要性，试验选取了5幅数据，对比了不同方法和步骤对结果的影响，如表 3所示。在表格中，“SRTM1_A”表示利用原始ASTER GDEM v2填充SRTM1空洞的结果，而“TIN_DSF”表示不经过精度校正，采用TIN差分曲面算法的结果，“MLR_DSF”和“ANN_DSF”分别表示由MLR模型和神经网络方法校正后，经TIN差分曲面算法填充的结果。从表 3可以看出，相比于SRTM1_A、TIN_DSF，MLR_DSF和ANN_DSF的RMSE值更小，说明精度校正能有效提高数据精度。不论是空洞还是整体的定量评价，经过校正的填充数据都能取得更优的结果。而对比ANN_DSF和MLR_DSF可以发现，本文方法在定量评价上优于MLR方法，表明其具有较强的空间预测能力，能更好拟合数据之间的关系，与表 1-2中点面校正模型的验证结果趋势相一致。

2.2 试验结果与分析

本文随机选取了不同空洞分布特征、不同空洞率的13幅数据进行试验，在图 5中给出了试验数据的空洞率和空洞区域内ICESat GLAS数据的分布。ICESat GLAS点的分布直接关系着评价的可靠性。图 6展示了SRTM3 v4.1、ASTER GDEM v2、SRTM1_A与本文方法融合结果的RMSE对比情况。其中，非空洞区域的评价反映不同DEM数据集的精度水平，空洞区域的评价体现各个方法的优劣。

从图 5可明显看出，空洞率和空洞区域内ICESat GLAS数目的趋势相同，说明ICESat GLAS在参考数据覆盖范围内沿轨道分布较为均匀。结合图 6可看出，在大多数情况下，融合结果的RMSE值都表现出明显的优势。在空洞区域内(图 6(a))，多数ASTER GDEM v2的垂直精度优于SRTM3 v4.1，而不及融合结果，可见本文方法对多数空洞区域的校正有效。对于“N28E101”和“N29E102”，SRTM3 v4.1的精度表现异常，原因在于地形陡峭，该区域的填充精度不准确，也进一步验证了其填充结果精度不稳定的现象^[28]。从图 6(b)中可以看出，大多数数据非空洞区域的垂直精度明显优于空洞区域。而“N29E87”和“N29E88”数据空洞范围非常小，因此验证样本在空洞区域数量较少，评价结果不具有代表性。总体而言，融合结果(图 6(c))的精度都优于其他DEM数据，在得到空间无缝DEM的同时，还能保证其具有稳定的垂直精度。

从视觉效果上来看，图 7展示了“N28E101”不同区域的放大图，该完整数据块的覆盖范围为28°N-29°N，101°E-102°E，空洞率达16.96%。图中，空洞范围的ASTER GDEM v2地形纹理较为清晰，为融合结果奠定了基础，但受噪声和异常值的影响。相比于ASTER GDEM v2和SRTM1_A，本文方法消除了噪声和异常值，且空洞填充的边界过渡也更加平滑，尤其在深色选框内效果显著。图 8展示了“N29E102”的区域彩色渲染图，该完整数据块覆盖范围为29°N-30°N，102°E-103°E，空洞率达10.04%。由图可知，该区域地形较崎岖且空洞范围大，从而影响地形的连续性。经过精度校正和填充，融合结果既能保留原始SRTM1的优势，又能有效地挖掘辅助数据的特征潜力，在目视判别上明显优于原始数据。

3 结论

本文结合高差拟合神经网络和TIN差分曲面方法，融合多源数据对SRTM1数据中的空洞进行修复。多源数据融合的方式既能保留SRTM1的细节信息，又能有效结合ASTER GDEM v2的无缝特性和ICESat GLAS的高精度特点，生成高质量的空间无缝DEM。通过试验发现，与不同方法相比，本文方法的结果具有最小的RMSE值。从视觉效果上来看，融合结果在填充区域具有更加自然的地形纹理，在填充边界过渡也更加平滑。目前，此方法对空洞中较大范围的噪声或异常值的处理存在一定的局限性。下一步可对辅助数据进行扩充，融入更多的高精度测高数据，如Cryosat2及即将发射的ICESat2数据，对方法进行改善，并考虑利用无参考修复的方法对残留异常值进行进一步处理。