基于ICESat-2数据及TanDEM-X DEM的林下地形反演

张 晨,朱建军,付海强

(中南大学 地球科学与信息物理学院，湖南 长沙 410083)

数字高程模型(Digital Elevation Model, DEM)作为最重要的基础地理信息产品之一[1]，已经在基础工程建设、水资源管理、地质学、冰川研究、气象以及军事制导等研究领域发挥着重要的作用[2]。随着对地观测技术的发展，目前获取的全球DEM向高精度、高分辨率、多样性的方向发展。然而受现有对地观测技术穿透性的限制，获取的DEM数据中包含植被、建筑物等高度信息。尤其在植被覆盖区域，传统光学遥感及微波遥感手段获取的是植被冠层的高程信息，并非准确的林下地表高程。为此，如何获得准确的植被高，将其从上述DEM产品中移除，对于获取大范围、高精度地表高程信息至关重要。

现阶段，主要的DEM产品有SRTM DEM[3]、ASTER GDEM[4]、AW3D30[5]及TanDEM-X DEM[6]。其中SRTM及TanDEM-X DEM分别通过C、X 波段InSAR (Synthetic Aperture Radar Interferometry)技术获取，ASTER GDEM及AW3D30产品通过光学遥感手段获取。在森林植被覆盖区，这些DEM受全部或部分植被高度信号的干扰，无法反映林下地形真实情况。因此，要想得到林下地形，需从上述DEM产品中扣除植被高信号。目前，现有研究主要是对SRTM DEM中包含的植被高信号进行扣除。2007年，Wilson等通过实地调查植被高结合Landsat-7遥感影像估算C波段SAR信号的穿透深度为植被高的50%。而后，从SRTM剔除相应的植被高[7]。2012年，Gallant等利用Landsat TM数据进行地块分析，通过最小二乘估算出地块边缘的植被高偏移量，并将其从SRTM DEM中扣除。但是，以上两种方法仅限于小范围，无法满足大范围林下地形制图的需求[8]。2013年，Baugh等人利用ICESat-1GLAS数据测量结果发现：C波段SAR信号在亚马逊热带雨林区的穿透深度为植被高的40%，随后，结合全球植被高产品，从SRTM中扣除植被高度。尽管该方法适用于大范围林下地形制图，但对C波段SAR信号穿透深度的计算方法过于简略，没有充分考虑穿透深度在空间上的显著差异[9]。针对该问题，2016年，O'Loughlin等利用ICESat-1GLAS、MODIS数据以及全球植被高产品，建立SRTM DEM中C波段SAR信号穿透深度与MODIS植被连续场(Vegetation Continuous Field, VCF)的关联，之后，利用VCF预测穿透深度的空间变化，最终实现SRTM DEM中植被高信号扣除[2]。

上述研究工作主要利用Simard等[10]发布的全球植被高产品。但是，该全球植被高产品所用数据为2007年之前获取，受植被高时空变化影响，该全球植被高产品无法用于TanDEM-X DEM(2018年发布)等新一代全球DEM中植被高信号的扣除。2019年5月NASA发布了ICESat-2/ATLAS(Advanced Topographic Laser Altimeter System)数据，其采用先进的光子点云数据采集技术，与ICESat-1/GLAS数据相比，ICESat-2/ATLAS具有更高的空间分辨率，对坡度具有更强的抵抗力，更易于测量山区地形以及植被高[12]。因此，如何利用ICESat-2/ATLAS数据获取大范围植被高并将其从TanDEM-X DEM中扣除，对提升DEM质量具有重要意义。基于此，本文提出一种基于ICESat-2/ATLAS数据和TanDEM-X DEM的林下地形生成方法，该方法主要包含三个步骤：①融合ICESat-2/ATLAS、 MODIS、温度、降雨及高程辅助数据并结合随机森林机器学习法获取植被高；②基于得到的植被高数据、TanDEM-X DEM及VCF数据，利用随机森林方法对X波段SAR数据的穿透深度进行建模及预测；③从TanDEM-X DEM扣除有效的植被高，得到大范围、高精确的林下地形。最后，选用我国东北区域作为实验区对所提算法的有效性进行验证。

1 扣除TanDEM-X DEM中植被信号的林下地形反演

1.1 基于ICESat-2 ATLAS植被高度反演

NASA(National Aeronautics and Space Administration)于2018年9月15日发射ICESat-2卫星，用于测量冰盖、冰川、冰架、海冰以及植被高度。ICESat-2搭载先进的地形激光高度计系统(ATLAS)传感器，通过采集光子点云数据对地表进行测量，具有小光斑、多波束、高采样率的特性[11]。采集的数据具有较高的空间分辨率，能记录丰富的地形信息，对坡度影响具有较强的抵御能力，在林下地形测绘及植被冠层探测等方面具有巨大潜力。本文主要采用ICESat-2/ATLAS ATL08产品中的地形及植被高数据，该数据在沿轨方向采样间隔为100 m[12]。由于地面和植被有较低的反射率，返回信号大部分来源于强光波束，因此需要根据升降轨来判别强弱光，从而选取强光信号；其次，数据采集期间，会有大量的噪声光子存在，但是夜间数据存在的背景噪声明显小于白天，这里需要根据夜间标志(night_flag)对数据进行筛选、剔除；最后，需要根据云标志(cloud_flag_asr)以及植被标志(canopy_flag)对数据进行筛选，从而得到可靠数据。但是，基于ICESat-2/ATLAS数据提取的植被高度仍为空间离散点，需考虑如何借助其他辅助数据填补其空白区。

从空间大尺度角度来看，影响植被分布的因素主要有纬度变化、海陆分布以及海拔高度三种情况。纬度影响主要体现在不同纬度温度不同，进而影响植被高度分布，即植被的纬度地带性；海陆分布影响主要是从海洋向内陆水汽含量逐渐变少，从而影响降雨量的变化，继而影响植被高度分布，即植被的经度地带性；海拔高度影响主要表现在随着海拔的变化温度也随之发生改变，使得植被分布随着海拔的增加呈现垂直分布，即植被的垂直地带性[13-14]。基于上述考虑，本文选取全球降水测量(Global Precipitation Measurement，GPM)数据、陆地地表温度产品MOD11A2、植被连续场(Vegetation Continuous Fields，VCF)数据MOD44B以及高程数据TanDEM-X DEM几个辅助数据进行植被高的反演。主要原理为通过寻求植被高度与上述辅助数据之间的关联，利用得到的关联及辅助数据推演空间连续植被高度图。

所有辅助数据均被采样至1 km的分辨率，如果一个分辨单元内包含多个ICESat-2点，则计算这些点的均值，以此表征该像元内的植被高。由于ICESat-2 ATLAS数据的精度受坡度变化的影响较小，因此，根据地面坡度对数据进行筛选。此外，尽管ICESat-2相较于ICESat-1具备较高的沿轨方向分辨率，但作为一种“点”数据，其难达到“面”测量的需求。为解决这一问题，本文利用ICESat-2/ATLAS(ATL08)植被高数据及上述辅助数据，利用随机森林(Random Forest， RF)进行训练，寻求二者之间的关联，在此基础之上，以辅助数据为输入，利用RF对植被高度进行预测[11]。

1.2 扣除TanDEM-X DEM中植被高信号的林下地形反演

当得到植被高度后，需要将反演的植被高产品及辅助数据进行重采样，使其与TanDEM-X DEM具备相同的分辨率(90 m)。TanDEM-X DEM通过X波段InSAR技术获得，由于X波段SAR信号具有一定穿透能力，测得的高程并非位于冠层顶部，不能从TanDEM-X DEM中直接扣除植被高(Hveg)，需要考虑穿透深度的影响。进而TanDEM-X DEM中扣除的植被高度Hrm可表示为：

Hrm=a×Hveg.

(1)

式中：a为扣除后植被高与估算的植被高度的比例系数。随后，通过移除DEM中的植被高度得到林下地形的高度信息，可以表述为：

DEMbare-earth=DEM-Hrm.

(2)

联合式(1)和式(2)，有效植被占比可以表示为：

(3)

植被覆盖区，InSAR穿透深度的主要与植被高、植被密度、地形、含水量等因素有关。为了准确估计每一个像素的扣除植被比，本文首先依据ICESat-2 ATL08产品提供的“裸”地表高程信息，计算其对应的扣除植被比；然后，通过RF方法构建扣除植被比a与植被高(Hveg)、VCF、高程(TanDEM-X DEM)之间的关系，依据所得到的模型，估算其余像素的扣除植被占比；最后，计算每一个像素对应的有效植被高度，并将其从TanDEM-X DEM数据中去除，即可得到林下地形信息。具体流程如图1所示。

图1 扣除TanDEM-X DEM中植被高信号的林下地形反演流程图

2 实验结果与分析

为了验证算法有效性，本文选取我国东北区域作为实验区，主要包括黑龙江省，吉林省，辽宁省，内蒙古东部(呼伦贝尔市，兴安盟，通辽市和赤峰市)以及河北北部(承德市，秦皇岛市，唐山市以及张家口市)。该区域属于温带季风气候，夏季高温多雨、冬季寒冷干燥。植被覆盖面积广阔，以针叶林和针阔混交林为主。ICESat-2数据分布及DEM图2所示。

图2 实验区及ICESat-2数据分布与DEM

2.1 植被高反演结果

基于图1算法流程，本文通过对ICESat-2 ATLAS数据以及辅助数据的处理提取计算、激光光斑植被高建模，使用随机森林算法外推出没有ICESat-2 ATLAS数据区域的植被高，得到研究区域植被高分布图，如图3所示。其中，大部分植被高在13 m左右，高植被区域主要在吉林省的东部、辽宁省的南部以及黑龙江省伊春市。占据植被分布面积较多的大兴安岭区域植被高主要集中在13～20 m。

图3 平均植被高图

对上述反演结果，用TanDEM-X森林/非森林(TanDEM-X Forest/Non-Forest，TanDEM-X FNF)数据进行掩膜，最后采用ICESat-2 数据中预留的少部分数据进行最终的精度评定。植被高制图结果一个分辨单元内若存在多个验证数据点，则计算其平均值，用该平均值代表该像元的平均植被高。如图4所示，本文反演的植被高结果与验证数据具有较好的一致性。植被高结果均方根误差RMSE为2.96 m，相关系数R2为0.87。误差主要分布在±5 m之间。反演结果存在低矮植被区域高估，高植被区域略微低估的现象。造成上述误差的原因，主要包含采样尺度误差，模型误差以及地理定位误差。采样尺度误差主要体现在：为了使得数据保持相同分辨率，从而进行采样处理；模型误差主要是随机森林建模过程中的误差；地理定位误差主要体现在多幅遥感影像进行纠正时精度不一致导致的误差。

图4 植被高反演结果验证

本文还根据全球土地覆盖图GlobCover 2009对植被类型分类进行精度评定。研究区域主要有阔叶林，针叶落叶常绿林以及阔叶针叶混交林几种植被类型。由图5可知，几种林分反演结果与总体结果趋势相近。该区域，优势树种是针叶林，主要集中在大兴安岭地区，其反演精度也相对较高，均值14.21 m，RMSE为2.55 m，R2为0.86，相对精度为17.96%。其次是阔叶林(>5m)占据较大面积分布，主要集中在黑龙江省伊春市以及辽宁省的东南部，均值为14.12 m，RMSE为3.0 m，R2为0.85，相对精度为21.21%。最后是混交林，其反演精度与针叶林相当，植被高均值为12.51m，RMSE为2.53 m，R2为0.87，相对精度为22.19%。

图5 不同植被类型高度反演结果

2.2 林下地形计算结果与精度验证

利用随机森林建模计算出扣除植被的百分比，在植被高制图的基础上得到TanDEM-X DEM中包含的植被高信号，即需要扣除的植被高，如图6所示。将得到的结果从TanDEM-X DEM中移除，得到林下地形。为了验证算法的正确性，通过ICESat-2地面高程点(815 710个)，分别对改正前、后的DEM验证进行精度验证。其统计结果如图7所示。

图6 扣除的植被高图

图7 DEM结果验证统计结果

从图7(b)可以看出校正前DEM误差整体表现出负的偏态，说明TanDEM-X DEM中包含了明显的植被高信号。而经过改正后的DEM，由于移除了植被高的影响，其误差分布在0附近。为了进一步验证精度，分别计算了校正前、后的均方根误差(RMSE)和误差均值。校正后DEM的RMSE为9.14 m，相较于原始的DEM(RMSE=11.65 m)，其精度提高了约21.6%；且其误差大多集中在±20 m之间。误差均值也由-6.5 m改善至0.03 m。综上结果表明：本文所提的方法可以较好的校正现有DEM中包含的植被高度信号。

3 结束语

本文以东北区域为例，基于辅助数据和ICESat-2/ATLAS数据，采用随机森林建模方法，反演植被高。植被高结果RMSE为2.96 m。根据研究区域主要植被类型的不同，对不同类型植被进行精度评定，结果表明该区域的优势树种反演结果精度相对较好。其后，利用所得的植被高产品，采用随机森林建模，推算出研究区域的林下地形，RMSE为9.14 m，相比原DEM精度提高了21.2%。结果表明本文方法可较好的移除现有DEM中植被高信号，为快速生成大范围、高精度的DEM提供了参考。尽管星载激光雷达具备反演大范围植被高的能力，但分辨率过低，适用于大尺度分析，难以满足高分辨率林下地形制图的需求。如何进一步提高植被高反演结果的分辨率将是后期研究工作的重点。