地表覆盖分类对于地球系统科学、全球环境变化和可持续发展研究具有重要意义。传统研究大多采用光学等遥感图像作为土地覆盖分类的数据源。而探索地表覆盖分类新的数据源，如星载激光测高、SAR影像等，则可作为光学遥感图像的有益补充。

ICESat/GLAS是首台对地球进行连续全球观测的激光测高仪，通过以40 Hz的激光重复频率发射一束波长为1 064 nm激光束对地球表面进行探测，最终在沿卫星轨道地面上形成直径约65 m的激光足印和172 m的激光足印采样间隔^[1]。以ICESat卫星数据作为土地覆盖分类数据源的研究主要有四种：（1）基于回波波形参数分析，如Doung等人依据回波高斯分解后的波形参数，对荷兰实验区域内四类地表（高树/森林、城市、水和裸地/低植被）进行分类，分类结果总精度为73%^[2]。（2）基于机器学习的方法^[3-4]，如Ghosh在印度Uttarakhand实验地区，通过使用K均值，中心点聚类算法和模糊C均值法进行土地覆盖类型分类，分类准确率最高可达89.41%^[3]。（3）波形理论分类器^[5-6]，或结合机器学习的方法与波形理论分类器^[6]，如Li等人通过结合波形理论分类器和支持向量机分类器，对格陵兰区域的水体和海冰进行分类，分类结果的总体准确率达到98.21%^[6]。（4）联合遥感影像处理提高分类精度，如Liu等人通过联合使用ICESat卫星数据与Landsat多光谱图像，提高中国河南省实验区域内四类地表（耕地、林地、灌丛和水体/裸地）的分类精度，从单一数据源约85%分类总体准确率提升到91%^[7]。

ICESat于2010年退役，其第二代卫星ICESat-2已于2018年9月由NASA发射升空。ICESat-2/ATLAS能以10 kHz的激光重复率发射3对激光束（532 nm波长）对地球表面进行探测。每对波束都包含一个强波束和弱波束（能量比为4∶1），并在沿卫星轨道地面上形成直径约17 m的激光足印和0.7 m的激光足印采样间隔^[8-10]。ICESat-2采用了光子计数的新技术^[11-12]，在对地观测时具备更高的空间分辨率和更多的地表轮廓描述信息。因此，ICESat-2/ATLAS有望进一步开发星载激光测高仪在土地覆盖分类中的应用潜力。但由于ICESat-2卫星发射时间较晚，基于该卫星数据进行土地覆盖分类的应用研究较少。如Zhang等人根据ICESat-2卫星数据在冰雪覆盖区与裸地背景噪声的差异，提出了一种在高原山区中对雪域与裸地进行分类的方法，并在阿克赛钦和帕米尔高原附近的研究区域对该方法进行验证，验证结果表明该方法的总体准确率在93%以上^[13]；Liu等人基于随机森林和逻辑回归，在北加利福尼亚州和新墨西哥州西部森林区域，进行燃烧区和未燃烧区分类，随机森林和逻辑回归的分类精度分别达到83%和76%^[14]。

为此，文中提出了一种基于随机森林的ICESat-2卫星地表覆盖分类方法，创新性地采用ICESat-2卫星产品中光子数目、不同类型光子水平和垂直分布比例、信噪比、太阳条件、大气条件作为分类的输入，并在中国长三角地区开展了多类地表覆盖类型分类实验以进一步探索该类卫星激光测高在多类型地表覆盖分类领域的应用潜力。

从ATLAS发射出的小部分光子在完成首次穿越大气层后，经地球表面反射，再次穿越大气层返回进入ATLAS望远镜中，最终被仪器电子设备记录下来。在此过程中，来自532 nm波长的太阳光的背景光子可能会到达探测器，且其中的一些背景光子也会被ATLAS记录下来。ATLAS记录到达时间的任何光子称为光子事件，该记录过程与光子的来源无关。因此，ATLAS记录的光子事件的数量取决于地球表面的几何结构和反射率、太阳条件以及大气中的散射和衰减^[15]。

ICESat-2卫星官方通过划分一定间隔的沿轨距离，对所探测到的光子进行分析和处理。不同等级的产品具有不同的沿轨距离间隔，如ATL03为20 m^[16]，ATL08为100 m^[15]。文中所开展的地表覆盖分类的空间最小分辨率正是各级产品对应的单个沿轨距离间隔。通过分析分布于每个沿轨距离间隔内的探测光子在空间上的分布情况，可以识别出不同的地表覆盖。而每个沿轨距离间隔内有效探测光子（反射于地表的探测光子）的比例则将决定地表分类正确性。

ICESat-2卫星官方在产品发布前所进行的滤波工作，可以很大程度地抑制背景噪声光子的影响，但并非完全消除，特别当沿轨距离间隔内的有效探测光子数量过低时。为了更好地保证地表分类的准确性，文中依据了ICESat-2卫星产品中的信噪比标签，通过设定门限值对低信噪比（<3）的数据进行剔除。

同时，对每个沿轨距离间隔内的有效探测光子进行分析时，希望有效探测光子是均匀分布在整段沿轨距离间隔内。因此，文中对沿轨距离间隔内出现水平方向上超过20%有效探测光子缺失的数据进行了剔除。

如1.1节所介绍的，ICESat-2卫星数据对地表采用了光子级灵敏度的探测。同时，对不同类别的地表分类时是通过分析分布于每个沿轨距离间隔内的探测光子在空间上的分布情况而进行的。因此，在对不同地表进行分类时，关于分类器属性的选择，文中主要考虑以下几个方面。

首先，在每个沿轨距离间隔内，具有不同反射率的地表所返回的探测光子数是不同的。例如，雪和冰表面在532 nm波长的激光束处具有高反射率（典型的朗伯反射系数介于0.8和0.98之间），而树冠和地形表面的反射率要低得多（土壤通常约为0.3，植被约为0.1），预计从地球表面返回的光子数量是冰雪表面的1/3到1/9^{[9, 15]}。基于此，文中采用产品在每个沿轨距离间隔内所记录的光子数量作为分类器的属性之一。该属性可以直接通过ICESat-2卫星产品中的标签获得。

其次，对于每个沿轨距离间隔，不同地表的探测光子在垂直方向中的不同高度层内，具有不同的比例。例如，当ICESat-2卫星对植被区域进行探测时，由于植被冠层的遮挡，分布于地面（或水面）的探测光子比例将低于裸地区域（或水域）；以及ICESat-2卫星对于森林区域与灌木等低植被区域的冠层内部和顶部进行探测时，由于植被遮挡情况的不同（树冠的茂密程度不同），两者的冠层内部和顶部具有不同的探测光子比例。基于此，文中在每个沿轨距离间隔内，计算被分类为地面光子、被分类为冠层光子、以及被分类为冠层顶层光子各自相对于全部光子的比例作为分类器的属性。上述前三类光子的数量可以直接通过ICESat-2产品中的标签获得。该类属性是关于地表垂直方向上的描述信息，为传统遥感影像所不具有的信息，将有望辅助传统光学遥感图像分类以进一步提高地表分类精度。

另外，对于每个沿轨距离间隔，不同地表的探测光子在水平方向中具有不同的分布特性（或比例）。例如，由于植被冠层的遮挡，裸地（或水体）区域内被分类为地面光子在水平方向上的分布相比于植被区域的地面光子更为均匀。基于此，文中主要在每个沿轨距离间隔内，计算被分类为地面的光子，以及被分类为冠层或冠层顶部的光子在水平方向上的分布比例作为分类器的属性。以ICESat-2卫星的ATL08产品为例，产品对于每个沿轨距离间隔，在水平方向上均划分为五段等长子间隔，并对每段子间隔内的光子类别标识从而构建指示标签^{[9, 15]}。通过此类标签，可以计算地面的光子，和冠层或冠层顶部的光子在水平方向上的分布比例。

最后，考虑到信噪比、太阳条件以及大气条件对有效探测光子的影响^[15]，文中将这些影响因素纳入分类器属性集内。这些属性可以直接通过ICESat-2产品中的标签所获得。

综上所述，文中构建分类器的属性如表1所示。

在完成数据预处理和分类器的属性选择后，通过选择训练样本进行随机森林训练可得到分类器模型，而测试样本可以通过该模型实现样本类别的预测，整个过程如图1所示。

图  1  基于随机森林的ICESat-2卫星数据地表覆盖分类流程图

Figure 1.  Flowchart of classifying land cover using ICESat-2 data with random forest

文中的实验区域位于中国的江苏省、上海市和浙江省，实验数据为ICESat-2卫星的ATL08产品^{[9, 15]}。产品中包含的IGBP土地覆盖分类数据可作为参考的地表覆盖数据^[15]。该数据是通过2012年的MODIS数据生成的，其空间分辨率为500 m，涵盖了从森林、灌木、草场、水体以及城市等17个不同类别。考虑到卫星激光测高在地表纹理、反射率等信息上表征的不足从而导致相似地表覆盖分类性能中存在局限性的情况，以及结合相关学者关于高树/森林、城市、水和裸地/低植被四类地表覆盖分类的研究^[2]，文中首先将水体、森林地表覆盖作为两类地表类别，然后从植被与非植被地表区分角度，初步对低植被，城市和裸地三类地表覆盖划分为低植被和城市/裸地地表覆盖作为后两类地表类别，最后依据上述四类地表覆盖分类结果的混淆矩阵，进一步合并混淆度较大的地表类别以进一步分析分类器性能。依据IGBP 17类地表的混淆矩阵^[17]，可以估算文中采用的水体、森林、低植被以及城市/裸地四类地表的总体分类精度约为93%。ICESat-2官方已对ATL08产品数据与该数据作了最佳的地理位置匹配^[15]，可以直接作为该实验的验证数据。

依据1.2节中所描述的数据预处理，首先通过ATL08产品中每个沿轨距离间隔内（100 m）的探测光子的信噪比标签（SNR），对信噪比低于3的数据进行剔除。然后，通过ATL08产品中的每个沿轨距离间隔内（100 m）的探测光子在水平方向上5段等长子间隔（20 m）类别分布标签（subset_te_flag），对超过20%有效探测光子缺失的数据进行剔除。剔除后的数据分布如图2所示。最后，将土地覆盖分类数据（IGBP土地覆盖分类数据）重新分为水体、森林、低植被以及城市/裸地四类。四种类别地表的样本数量如表2所示。

图  2  实验数据的分布

Figure 2.  Distribution of experimental data

依据1.3节，完成了属性集的构建，然后，对表2中的强波束和弱波束数据分别按以下步骤进行处理：

（1）分别从各个类别的样本中随机抽取25%作为训练样本。

（2）对训练样本进行训练，获得随机森林分类器。

（3）以剩余的样本作为测试样本，通过步骤（2）获得到的随机森林分类器进行对样本类别进行预测。

（4）重复步骤（1）~（3）五次，得到五组随机分类结果。

训练过程中，随机森林的参数中的ntrees设置为500，mtry设置为4。

按照2.2节的数据处理步骤进行处理后，可以获得五组强波束和弱波束数据关于四类地表的分类结果，如图3所示。其中，第一组强波束和弱波束各自的混淆矩阵如图4所示。

图  3  强波束和弱波束数据关于四类地表的分类结果

Figure 3.  Classification results for four kinds of land cover using strong beam and weak beam altimetry data

图  4  强波束和弱波束测高仪数据关于四类地表的分类混淆矩阵

Figure 4.  Classification confusion matrix of four kinds of land cover using strong beam and weak beam altimetry data

从图3中可以看出，强波束和弱波束的数据在水体、森林、低植被以及城市/裸地四类地表的分类结果较好，总体分类精度均优于85%，Kappa系数均优于70%。从图4中可以发现，低植被的分类精度最佳，往后依次为水体、森林和城市/裸地。其中，城市/裸地类别的地表与低植被类被的地表混淆最大（强波束：76.21%，弱波束：82.06%），这可能是因为部分低植被区域的冠层覆盖率较低，信号（在这种情况下是树冠）光子可能被错误地识别为数据产品内的太阳背景噪声而被过滤掉^[15]，从而与裸地（或城市中空旷区域）混淆。此外，城市中稀疏的建筑与植被交错区域也将是潜在的混淆源。为了进一步提高分类器的性能，笔者等尝试将低植被类别的地表与城市/裸地类别的地表作为相同类别的地表，从而开展水体、森林、以及低植被/城市/裸地三类地表分类实验。该地表分类实验结果如图5所示。从图5中可以看出，分类器的性能得到进一步的提高，强波束和弱波束的总体精度均提高至90%以上，Kappa系数均提高至79%以上。

图  5  强波束和弱波束测高仪数据关于三类地表的分类结果

Figure 5.  Classification results for three kinds of land cover using strong beam and weak beam altimetry data

文中提出了一种基于随机森林的ICESat-2卫星地表覆盖分类方法，并通过在中国长三角地区开展的多类地表覆盖类型分类实验对该方法进行了验证。验证结果表明，在水体、森林、低植被以及城市/裸地四类地表中，强波束和弱波束的总体精度均优于85%，Kappa系数均优于70%；水体、森林以及低植被/城市/裸地三类地表中，强波束和弱波束的总体精度均提高至90%以上，Kappa系数均提高至79%以上。该方法总体上所具有的良好分类性能，特别是分类过程中所兼顾的地表垂直方向上的描述信息（探测光子在垂直方向中的不同高度层内的分布比例），将有望辅助传统光学遥感图像分类以进一步提高地表分类精度。更多地，该方法在进行四类地表分类中，城市/裸地该类别的地表分类精度表现性能较差，主要是因为与低植被类别的地表混淆较大。这是在未来的工作中需要进一步研究改进的方向。