ICESat_2机载试验点云滤波及植被高度反演_夏少波.pdf

1、1007-4619( 2014) 06-1199-09 Journal of emote Sensing 遥感学 报ICESat-2机载试验点云滤波及植被高度反演夏少波1， 2， 王成1， 习晓环1， 骆社周1， 曾鸿程31 中国科学院遥感与数字地球研究所 数字地球重点实验室 ， 北 京 100094;2 中国科学院大学 ， 北京 100049;3 University of Toronto， Toronto， Canada摘 要 : 新一代星载激光雷达卫星 ICESat-2 将采用多波束微脉冲光子计数技术 ， 并进行高程剖面式的对地观测 。由于该点云数据具有背景噪声大 、密度低并呈线状分布等

2、 特点 ， 传统的点云滤波算法并不适用 ， 研 究 新的点云滤波算法十分必要 。本文以 ICESat-2 的机载模拟器 MABEL 数据为例 ， 首先介绍了微脉冲光子计数激光雷达的基本原理和数据特点 ， 并针对高程剖面点云提出基于局部距离统计和最小二乘局部曲线拟合的点云滤波算法 ; 然后 ， 对美国加利福尼亚州 Sierras-Forest 地区 MABEL 试验中 532 nm 通道的光子点云进行滤波处理 ， 并利用识别的地面点插值得到 3 m 分辨率的线状 DEM， 进而估算了该区域美国云杉的平均树高 ; 最后 ， 对该滤波算法进行精度评价 ， 并分析了误差来源及其对 DEM 精度和树高反

3、演精度的影响 。结果表明 : ( 1) 该算法整体精度达 97 6%， 能有效剔除绝大部分噪声点且对地形起伏具有较强的自适应能力 ; ( 2) 误分噪声点影响了滤波过程中局部地形的拟合 ， 而滤波过程中的分类误差将降低 DEM 和树高反演的精度 。关键词 : ICESat-2， MABEL， 微脉冲光子计数 ， 剖面点云 ， 去噪与滤波 ， 植被高度中图分类号 : TP79 文献标志码 : A引用格式 : 夏少波 ， 王 成 ， 习 晓环 ， 骆社周 ， 曾鸿程 2014ICESat-2 机载试验点云滤波及植被高度反演 遥感学报 ， 1 8( 4) : 11991207Xia S B， Wan

4、g C， Xi X H， Luo S Z and Zeng H C 2014 Point cloud filtering and tree height estimation using air-borne experiment data of ICESat-2 Journal of emote Sensing， 1 8( 4) : 1199 1207 DOI: 10 11834/jrs20144029收稿日 期 : 2014-02-10; 修订日期 : 2014-05-29; 优先数字出版日期 : 2014-06-05基金项目 : 国家自然科学基金重大国际合作研究项目 ( 编号 : 4112

5、0114001) ; 中国科学院百人计划专项第一作者简介 : 夏少波 ( 1991 ) ， 男 ， 硕士研究生 ， 研究方向为激光雷达遥感 。E-mail: shaoboxia1023 gmail com通信作者简介 : 王成 ( 1975 ) ， 男 ， 研究员 ， 研究方向为激光雷达遥感 。E-mail: wangcheng radi ac cn1 引 言新一代星载激光雷达卫星 ICESat-2( Ice， Cloud，and land Elevation Satellite-2) 将于 2017 年发 射 ， 其搭 载的 ATLAS( Advanced Topographic Laser

6、 Altime-ter System) 测高系统采用了微脉冲多波束光子计数激光雷达技术 ( Evans， 2013) ， 这是该技术首次应用于星载平台 ( Abdalati 等 ， 2010; Anthony 等 ， 2010) 。为了论证系统设计和进行算法预研究 ， NASA 模仿ATLAS 系统研发了 MABEL( Multiple Altimeter BeamExperimental Lidar) 激光雷达系统并进行了高空机载试验 ( Brun 等 ， 2011) ， 对其获取的数据进行处理及 应用研究对 ICESat-2 有重要意义 。与现阶段常见的脉冲检测式激光雷达 GLAS( Ge

7、oscience Laser Altimeter System) 、LVIS ( Land，Vegetation， and Ice Sensor) 和 ATM( Airborne Topo-graphic Mapper) 波形或点云类似 ， MABEL 记录光子返回处的 3 维坐标 ， 但其成像原理和数据特点与脉冲检测式激光雷达有所不同 ( Krainak 等 ， 2010) 。首先 ， MABEL 点云数据背景噪声高 ， 这是由微脉冲光子计数激光雷达的工作模式造成的 ， 该系统发射和检测的信号均为弱信号 ， 难以区分返回脉冲信号 、大气散射噪声 、太阳辐射噪声和仪器自身噪声( Harding

8、， 2009) ， 因此点云去噪算法在 MABEL 数据处理和应用中十分重要 。其次 ， MABEL 与传统的推扫式或扫描式激光雷达不同 ， 其数据沿飞行1200 Journal of emote Sensing 遥感 学 报 2014， 18( 6)轨迹 呈 窄带状分布 ， 属于高程剖面点云 ( McGill等 ， 2013) ， 且点云密度较低 ， 因此需要研究合适的点云滤波算法 。当前针对剖面式光子点云的去噪和滤波算法可分为两类 : 一是将剖面点云栅格化为 2 维影像 ， 每个像元包含了点云密度信息 ， 然后采用影像处理技术 ， 如 Canny 算子检测边界 ( Magruder 等 ，

9、 2012) 、中值滤波和轮廓检测 ( Awadallah 等 ， 2013) 等 ， 剔除噪声点并探测地面点 ; 二是计算点云剖面中每一个点的局部统计量 ， 利用统计量的分布特征 ( 直方图 ) 设置统计量的全局阈值并对点云进行逐点分类 ， 分为噪声点 、地面点和非地面点 ( Herzfeld 等 ， 2013;Brunt 等 ， 2014) 。前者需要对点云栅格化 ， 会造成信息损失 ， 降低滤波的精度 ; 后者需要针对统计量设置分类阈值 ， 而阈值的选择受到坡度 、地表类型 、点云密度和统计范围等因素的影响 ， 算法自适应性有待提高 。该类方法能保留更多的点云信息 ， 且在滤波分类前首先

10、剔除了大部分背景噪声 ，因此受噪声点的影响较小 ， 其在统计参量 、统计半径以及自适应阈值的选择方面均需要进一步研究 。目前已有的背景噪声剔除算法可以识别 90%的噪声点 ( Magruder 等 ， 2012) ， 但残留的噪声点比例相比非光子点云依然很高 。本文主要研究低背景噪声 ( 原始点云噪声点比例低于 20%， 或已进行过粗去噪处理 ) 的 MABEL 点云数据的滤波分类算法 ， 包括剖面点云中噪声点剔除和地面点与非地面点的准确分类 ; 研究区域位于美国加利福尼亚州 ( 加州 ) Sierras-Fores 林区 ， 在对 MABEL 数据滤波分类结果的基础上插值得到 DEM， 进而

11、估算平均树高 。2 MABEL 原理 与 试验21 微脉冲光子计数激光雷达基本原理按照 激 光雷达发射和接收信号的工作模式 ， 可将激光雷达系统分为脉冲模拟检测式和微脉冲光子计数式 ( Harding， 2009) 。前者发射的脉冲信号能量较高 、持续时间较长 ( 48 ns) ， 激光接收器累积一定时间间隔内返回的所有光子的能量 ， 并对累积能量的时间序列进行模数转换 ， 然后存储波形信息 。系统根据能量阈值对波形信息进行峰值提取及坐标计算 ， 从而产生离散点云记录 。由于这种模式需要高能量脉冲 ， 并通过在短间隔时间内累积大量光子能量的方式来保持接收信号较高的信噪比 ， 但却限制了仪器的发

12、射频率 ， 也加快了仪器的损耗速率 。相比而言 ， 微脉冲光子计数激光雷达发射的脉冲能量较低 、频率高 、持续时间短 ( 1 ns) ， 因而得名“微脉冲 ”。在记录返回信号时 ， 往往只需记录返回的数以万计个光子中的一个或几个的坐标 ， 无需光子数量的累积 ， 也无需提取波峰 ， 即可记录地表信息 ， 降低了波形分解及离散点云记录的不确定性 。但由于该方法检测的是弱信号 ， 故受噪声影响更大 ， 尤其是太阳背景噪声 。22 MABEL 设计与机载试验ICESat-2 上搭载的 ATLAS 系统将发射的激光分为 6 束 3 列 ， 并排列为面阵 ， 地面列间距约为 33 km， 且同一光束的激

13、光脚点沿轨间距约 0 7 m( 频率为 10 kHz) ， 脚点直径约为 10 m。MABEL 的激光器配置模仿了 ATLAS( 图 1) ， 区别在于 ， MA-BEL 工作波段为 1064 nm( 近红外 ) 和 532 nm( 绿色 ) ， 而 ATLAS 只有 532 nm。图 1 MABEL 多波束通道配置 ( McGill 等 ， 2013)Fig1 Multi-beam configuration of MABEL ( McGill， et al ， 2013)夏少 波 等 : ICESat-2 机载试验点云滤波及植被高度反演 1201MABEL 搭载 在 NASA 的 E-2

14、飞机上 ， 为接近ICESat-2 的实际工作情况 ， 其飞行高度约为 20 km，尽可能避免了大气干扰 ; 发射频率可在 525 kHz范围内调整 ， 当以 10 kHz( ATLAS 频率 ) 工作时 ， 脉冲宽度约为 2 ns， 地面脚印直径约为 2 m。MABEL从 2010 年开始飞行试验 ， 已经获取了覆盖森林 、农场 、戈壁 、沙漠和冰川等区域的数据 ， 产品格式为HDF5， 包含了每个光子的地理坐标 。表 1 列出了MABEL 及机载试验的相关参数 。表 1 MABEL 试验参 数Table 1 Parameters of MABEL experiment参数 指 标 值飞行高

15、度 约 20 km波段 532 nm/1064 nm激光光束数 215 个 ( 选择使用 )工作频率 525 kHz发射脉冲宽度 2 ns( 10 kHz)脚印直径 约 2 m记录 数 据 光子点云3 MABEL 点云 滤 波算法31 MABEL 数据 特 点分析与常见的点云滤波处理相比 ， MABEL 数据滤波的特殊性表现在 3 个方面 。( 1) 噪声大且无特定模式 ， 与观测环境密切相关 ， 而且难以预测 。背景噪声严重干扰点云分类 ， 因此去噪是点云滤波算法的重要组成部分 ， 简单的去噪算法 ( 如通过设置高程阈值等 ) 无法满足 MABEL 数据处理的需要 。( 2) 剖面式点云 ，

16、 其沿飞行轨迹方向呈窄带状分布 。现有的机载点云滤波算法如不规则三角网滤波 ( Axels-son， 2000) 、移动曲面滤波 ( 苏伟 等 ， 2009) 等需要构网或曲面拟合 ， 因此并不适用于接近线状分布的点云数据处理 。( 3) 点云密度低 ， 如本文采用试验数据的点云密度约 0 87 pts/m2， ICESat-2 设计 的 沿轨道采样密度也不超过 1 5 pts/m2， 给滤波算法中的种 子 选 取以及坡度滤波 ( Vosselman， 2000) 中角度 、边长 、高差等阈值选择带来困难 。因此 ， 目前常用的点云滤波算法不能直接应用于 MABEL 及 ICESat-2 的点

17、云数据处理中 。本文滤波算法框架包含了点云精去噪以及地面点与非地面点分类 ， 并主要采用了自适应阈值方法 ， 算法流程如图 2。首先将点云投影到局部坐标系 ， 并将点云的平面坐标转换为沿轨距离 ， 从而生成 2 维高程剖面点云 ; 然后利用点云精去噪算法进一步剔除数据中残留的噪声点 ， 并采用 K-D Tree 结构建立点云的分段索引 ， 通过局部最小二乘曲线拟合方法对点云进行滤波分类 ; 最后基于滤波分类结果 ， 采用曲线拟合方法插值得到 DEM。图 2 MABEL 剖面点云滤波算法流程Fig2 Flowchart of MABEL profile filtering algorithm32

18、 基于局部距离统计的点云精去噪光子计 数激光雷达输出的点云包含了地物反射 、大气散射 、太阳背景噪声和仪器噪声 。由于 MA-BEL 数据的大气散射和仪器噪声的能量比地物反射与太阳背景低很多 ， 一般可以忽略 ( Kwok 等 ，2014) ， 因此 ， 光子点云去噪研究的主要对象是太阳背景噪声 。太阳背景噪声与非噪声点在空间分布特征上区别较大 ， 噪声点的局部密度更低 ， 其与临近点的距离比非噪声点与临近点的距离更大 ， 一般相差 12 个数量级 。本文考虑了噪声点与非噪声点空间分布特征的差异 ， 并基于室内场景 2 5D 点云去噪方法1202 Journal of emote Sensin

19、g 遥感 学 报 2014， 18( 6)( usu， 2009) ， 设计了适合光子剖面点云的精去噪算法 。算 法 主要包含 3 个步骤 :步骤 1 计算点云中每个点到周围最临近的 k个点的总距离 disti。disti=j = kj =1( xi xj)2+ ( hi hj)槡2( 1)式中 ， x 表示沿轨道距离 ， h 表 示 点云高程 。步骤 2 统计局部总距离的频数直方图 。步骤 3 设置阈值剔除噪声点 。假设无噪声点点云的 disti的频 数 直方图呈高斯分布 ， 其均值与方差通过点云局部统计进行估算 。若 disti大于 其 整体均值与 t 倍标准差之和 ， 则 i 点为噪声点

20、 。将累计距离直方图近似为高斯分布的前提假设是 “区域内非噪声点云分布较为均匀 ”。若非噪声点云近似均匀分布 ， 而样本点数目较大 ， 则点云中每个点到其周围一定数目点的累计距离将十分接近 ， 距离和远大于或小于全局平均距离的点很可能是孤立点 。点云空间分布不仅与激光器采样方式有关 ， 还与地物类别有密切联系 。以植被和建筑为例 ， 植被点云不仅仅分布在冠层顶部 ， 还分布在冠层内部以及林下 ， 而建筑物的点云只分布在表面 ， 并随建筑物的几何外形而变化 。因此 ， 不同地物类别的点云空间分布模式 ( 点间距和点密度等 )不同 ， 混合场景的距离分布频数直方图可以视为多个高斯分布的叠加 。本文

21、的研究对象为森林 ， 且沿轨采样比较均匀 ， 故可将非噪声点云的累计距离直方图近似为高斯分布 。k 值是计算距离和的核心参数 ， 即参与计算的临近点个数 ， 实质上是一个平滑因子 ， 其作用是平滑同类别点云 ( 非噪声点 ) 内部的空间异质性 。若 k值过小 ， 则平滑作用不明显 ， 其频数直方图或可表现为多个高斯分布的叠加 ; 若 k 值过大 ， 则平滑了噪声点云与非噪声点云之间的空间分布差异 ， 使点云去噪难以进行 。在对同一地物类型 、不同区段的剖面点云进行试验和分析后 ， 选择 k =50 作为经验值 ，但对于非植被区域或点云密度与本文试验数据相差较大的区域 ， k 值需要重新调整 。

22、由于存在较多噪声点 ， 其累计距离分布影响了高斯参数的计算 ， 故本文对频数直方图进行高斯拟合时 ， 频数峰值对应的距离作为近似均值 ， 并将峰值位置到最小累计距离频数位置的距离作为近似标准差 ， 且 t =20。33 基于最小二乘局部曲线拟合的点云滤波本文采用局部曲线拟合的方法对非噪声点进行了地面点与非地面点的分类 。该方法与基于局部统计直方图的 MABEL 滤波算法 ( Herzfeld 等 ，2013; Brunt 等 ， 2014) 不同 ， 有更明确的几何意义 。其本质是将曲面滤波中的曲面拟合过程褪化为曲线拟合 ， 分为 4 个步骤 :步骤 1 查找局部最低点 。利用 K-D Tre

23、e 建立点云索引 ， 在一定大小的窗口内查找高程最低点 。若存在较长数据段中地面点缺失 ， 则修改查找窗口重新查找 。步骤 2 利用局部最低点拟合局部地形 。对步骤 1 中的每个局部最低点进行临近搜索 ， 联合与其沿轨距离最近的 4 个局部最低点 ， 采用最小二乘算法求解局部二阶曲线 。h = ax2+ bx + c ( 2)式中 ， x 表 示 沿轨距离 ， h 表示高程 ， 参数 a、b 和 c 为待求的拟合参数 。步骤 3 计算自适应阈值 。地面点与非地面点通过高差阈值进行区分 。为了使算法更具普适性 ，采用自适应阈值 。阈值的选择参考了机载点云移动曲面滤波算法 ( 苏伟 等 ， 200

24、9) :Threshold = s ( Hmax Hmin) ( 3)式 中 ， Hmax表示局部窗口内点云高程最高值 ， Hmin表示局 部 窗口内点云高程最低值 ， s 可依据地形起伏设置 ， 本文涉及的实验中 s =01。步骤 4 点云分类 。利用步骤 2 中局部地形曲线求得沿轨各点的理论地形高程 ， 激光脚点高程减去理论地形高程得到两者高差 ， 若某点高差大于步骤 3 中计算的自适应阈值 ， 则该点判定为非地面点 ，反之判定为地面点 。4 实验 与 分析41 实验过程与结果实验 区域位于美国加利福 尼亚州 Sierras-Forest地区 ( 图 3( a) ) ， 该 区域主要植被为

25、美国云杉 。采用2012 年 9 月飞行试验的 532 nm 波段第 11 通道的观测数据进行实验 ( 图 3( b) ) 。原始实验数据如图 3( c) 所示 ， 其横坐标是时间 ( s) ， 纵坐标表示高程 ( m) 。实验数据中存在部分噪声点 。由于该段数据跨度长 ，点云密度变化大 ， 故本文选取 512 s 数据段进行实验 ， 其不仅包含较为平缓的地形 ， 还存在一段长度为2 00 m、坡度约为 11的斜坡 ， 具有一定代表性 。首先将光子点云进行坐标投影变换并沿地面飞行轨迹方向排列 ， 进而生成高程剖面点云 。如图夏少 波 等 : ICESat-2 机载试验点云滤波及植被高度反演 1

26、2034所示 ， 横坐标表示沿轨距离 ， 纵坐标表示点云高程 。试 验 数据段总点数为 2001 个 ， 沿轨距离约为2310 m。图 3 实验区域与点云剖面Fig3 Location of the study area and the profile of points cloud图 4 实验样本原始点云Fig4 Original data of the study area采用 31 节中的去噪算法对点云进行精去噪 ， k值设为 50， 即累计该点与临近 50 个点的距离 。disti的分 布 如图 5( a) 所示 ， 并统计其频数直方图如图 5( b) 所示 ， 其大多数距离和在 50

27、 m 以内 ， 图 5( b) 中的子图表示距离和大于 50 m( 绿色矩形内 ) 的局部放大 。图 5( b) 中 ， 局部距离和在 050 m 的频数近似呈高斯分布 。频数直方图的峰值对应的距离为 225 m， 而最小累计距离为 14 8 m， 即近似均值为22 5 m， 近似标准差为 7 7， 当 t =2 0， 则距离和大于 379 m 的点归为噪声点 ， 并予以剔除 。精去噪前后的点云如图 6 所示 ， 高空及地表下的噪声点基本被剔除 ， 其中判别为噪声点的总数为 82 个 ( 红色点云 ) ， 剩余非噪声点共 1919 个 ( 绿色点云 ) 。图 5 点云精去噪阈值选择Fig 5

28、Threshold selection for points cloud denoising图 6 精去噪算法识别的噪声点 ( 红 ) 与非 噪 声点 ( 绿 )Fig6 Noise points ( red) and no-noise points( green)detected by denoising algorithm接着将局部最低点查找窗口大小设为 30 m， 以确 保最低点不落在植被 ， 尤其是冠幅较大的植被上 。图 7( a) 中红色点云即选择的局部最低点 ， 总点数为 77 个 。进而搜索每一个局部最低点的 2 维临近点 4 个 ， 利用最小二乘方法拟合局部二阶曲线 ， 并依据

29、 32 节方法计算局部阈值 ， 然后依次对所有激光脚点进行分类 ， 结果如图 7 ( b) 所示 ， 共计 760 个1204 Journal of emote Sensing 遥感 学 报 2014， 18( 6)点被判定为地面点 ， 1159 个点被判定为植被点 。将线 状 DEM 分辨率设为 3 m， 利用临近待求格网点的 20 个已分类地面点进行二次曲线拟合 ， 从而插值得到格网点高程 ， 如图 7( c) 所示 。最后 ， 假设树冠冠幅为 15 m， 沿距离方向求出 15 m 间隔区间内点云的高程最大值作为树冠点 ， 如图 7( d) 所示 。将树冠点高程值减去 DEM 对应位置高程

30、值得到树高 ， 进而计算出实验区域的平均树高为 342 m。图 7 地面与植被分类过程及结果Fig7 Denoising and filtering results42 分析 与 讨论对点云滤波效果的分析方法 ， 可以结合高精度高分辨率的 DEM， 而 MABEL 试验区点云的空间分辨率约为 1 m， 且其观测精度能达到厘米级 ( Brunt等 ， 2014) ， 因而难以获取满足如此高分辨率和高精度的地形数据 。同时 ， MABEL 点云不仅采样频率低 ， 脚印大 ( 1 m) ， 而且混入很多噪声点 ， 因此 ， 点云自身即含有较大的不确定性 。即使能结合真实的参照标准 ， 如高精度地形

31、、实地测绘等 ， 临近地面点与非地面点的分类也存在问题 。总之 ， 结合 MA-BEL 数据特点和现有条件 ， 绝对意义的分类精度评定存在困难 。本文将分类结果和人工判读结果进行比较 ， 主要目的是辅助分析和评价算法的误差来源和影响程度 。需要说明的是 ， 人工判读光子点云时 ， 那些与地面或者植被点云接近的激光点依旧难以准确划分 ， 并且也不能排除冠层内包含有噪声点的可能性 ， 但这是目前最可行有效的算法精度分析途径 。统计得到的混淆矩阵如表 2。表 2 点云分类结果的混淆矩阵Table 2 Confusion matrix of the classification result /个地面

32、 植 被 噪声 累计地面 726 14 0 740植被 30 1145 0 1175噪声 4 0 82 86累计 760 1159 82 2001分析 混 淆矩阵可知 ， 算法的整体精度为 976%。其中 ， 植被点与地面点错分的数目最多 ， 不但存在地面点被错分为植被点 ， 而且也存在植被点被误认为地面点的情况 。其次 ， 精去噪误差主要来自地面点与噪声点分类 。下面结合误差来源对滤波算法与植被高度反演进行分析与讨论 。421 点云精去噪结果分析表 2 显示点云精去噪算法识别了 82 个噪声点 ，同时 4 个噪声点误判为非噪声点 ， 但是没有任何非噪声点误判为噪声点 。错分主要有两个原因 :

33、( 1) MABEL光子点云的噪声散布在传感器下方的空间内 ， 而不局限在地物附近 ( Herzfeld 等 ， 2013) ， 因此 ， 空中或者远低于地表的噪声点的局部累计距离与地物附近的噪声点的局部累计距离之间的量级差异较大 ， 空中散布的噪声点其局部累计距离可达数百米 ， 而近地表的局部累计距离往往低于 100 m;( 2) 由于 MABEL 点云较为稀疏 ， 接近地表但低于实际地面的噪声点其局部累计距离与地表以上植被点的局部累计距离接近 ， 即两者的空间邻域特征近似 。当确定 t 值后 ， 由于局部的阈值是唯一的 ， 其不能同时分辨出非噪声点 、较高空噪声点及与非噪声点接近的噪声点

34、。将地物附近的噪声点与地物点进行分离 ， 对于接近地表的噪声点 ， 需要有真实地面高程作为分类依据 ， 而真实地面高往往是未知的 ; 而与地物接近的地表之上的噪声点 ， 由于点云稀疏 ， 人工尚且无法准确判读 。因此 ， 这类噪声的剔除需要更深入的研究 。夏少 波 等 : ICESat-2 机载试验点云滤波及植被高度反演 1205422 点云滤波结果分析结合 表 2 与 图 7 发现剖面点云分类的效果很好 ， 基本区分了地面点与植被点 。部分数据段如850950 m 段 ， 植被较密而导致地面点缺失 ， 但是由于滤波窗口选择合适 ， 该类点云并没有被错分 。在沿轨 8001000 m 范围内

35、， 地形坡度约为 11， 而滤波算法生成的地形曲线较好地拟合了该长坡 ， 说明本文算法具有一定的自适应性 。混淆矩阵表明植被点和地面点存在较多相互错分的点 ( 共 44 个 ) ， 结合图 7( b) 可以发现这些误判点主要分布在地形点与植被点接近的位置 。这类点的错分主要是由地形拟合误差和滤波阈值设置不当引起的 ， 而地形拟合与阈值计算均与局部种子点的选取密切相关 。事实上 ， 原始点云进行去噪处理后 ， 依然残留了 4 个噪声点 ， 若这些噪声点被选为地面种子点 ， 即有可能造成分类错误 。以图 8 例 ，沿轨约 800 m 处的点 ( 黑色箭头指向 ) 人工判定其为地表下的噪声点 ， 精

36、去噪算法将其判别为非噪声点 ， 因此在寻找局部最低点时该点被误判成地面点 ， 若采用该点进行局部曲线拟合 ， 将导致拟合曲线比真实地形曲线低 ( 图 9( a) ) ， 且计算的自适应阈 值也 过大 ， 故真实地面点 ( 红色矩形框内的点 ) 被错判为植被点 。由此可见低于地表的噪声点直接影响了点云滤波分类的结果 。图 8 图 7( b) 局部 放 大 ( 黑色箭头指向噪声点 ，红色矩形框包含了误分的地面点云 )Fig8 Partial enlarged drawing of Fig7( b)( The black arrow points to the noisy return， the r

37、eturns withinthe rectangle indicate the misclassified ground points)如果需 要在拟合中剔除地表以下的噪声点 ， 那么可以采用类似迭代曲线拟合的方法 ( 王峰 等 ，2011) 。仍以图 ( 8) 中黑色箭头所指向的噪声点为例 ， 其附近 5 个初选地面点的二阶曲线拟合结果如图 9( a) 所示 。图中红色点即噪声点 ， 导致曲线向噪声点弯曲 ， 且噪声点的拟合残差达到 3 5 m( 图 9( b) ) ， 若 采用 一定准则将该点剔除 ， 然后进一步对曲线进行拟合 ， 或可得到合理的地形曲线 。然而 ，由于 MABEL 点云密

38、度较低 ， 同时为了保证地面点选取的精度 ， 种子地面点的分布较为稀疏 ， 因此若采用稳健估计的迭代方法 ， 其阈值选取较难 ， 且由于数据的抽稀 ， 可能会影响地形曲线的精度 ， 进而影响分类结果 。因此 ， 这类方法仍需进一步研究 。图 9 含噪声点的局部曲线拟合及其残差分布Fig9 Error distribution by local curve fitting with noise423 滤波 误 差对 DEM 及树高估算的影响由于存在错分点 ， 尤其是低于地表的噪声点 ，因此 ， 由地面点插值生产线状 DEM 时 ， 必然引入误差 。实验发现联合临近多点的 DEM 曲线插值可以削弱

39、局部点云错分的影响 。如图 10， 绿色点为原始点云 ， 由于噪声点存在 ， 导致红色矩形框内的地面点错分为植被点 ( 图 ( 8) ) 。图 10 为通过联合局部20 个已分类的地面点进行二阶曲线拟合 ， 并内插得到 DEM 格点序列 ( 灰色方块 ) 。噪声点 ( 黑色箭头指向点 ) 出现在 DEM 格点序列下方 ， 而错分为植被点的地面点落在了 DEM 格点之上 ， 这表明局部点云错分可通过一定数目的临近地面点进行纠正 ， 但只适用于错分地面点较少的情况 ， 并且图 10 中红框内DEM 格点依旧比实际地表 ( 红框内 5 个绿色点 ) 低 。由于缺乏高精度地形数据进行对比 ， 因此 ，

40、 插值生1206 Journal of emote Sensing 遥感 学 报 2014， 18( 6)成线 状 DEM 的误差无法进行定量分析 。由于目前没有实测的树高 ， 因此树高反演误差无法定量评价 ， 本文将其作为应用实例进行分析 。与 ICESat-1 通过分析大脚印波形参数 ( 董立新 ，2008) 来反演植被高度不同 ， MABEL 以及未来的ICESat-2 将采用类似于机载激光雷达反演树高的方法 ( 庞勇 等 ， 2006) ， 通过 DSM 和 DEM 的差值反演植被高度 。即将一定间隔内的最高点作为树顶点 ，树顶点高程减去对应的 DEM 高程得到了平均树高 。滤波误差影

41、响 DEM 生产误差 ， 进而影响反演树高的精度 。图 10 图 7( c) 局部 放 大 ( DEM 格点 ( 灰色方块 ) 与去噪点云叠加 ，黑色箭头指向残留噪声 ， 红色矩形框则包含误分的植被点云 )Fig10 Partial enlarged drawing of Fig7( c) ( The grey blocksare DEM grids， the black arrow points to the noisy return，the red rectangle indicates the wrongly labeled as ground points)5 结 论本文 提 出的剖面

42、点云滤波算法简单快速 ， 具有较强的坡度自适应性 ， 且局部分段拟合过程有利于长轨迹数据的分段并行处理 。由于光子点云的噪声很难完全剔除 ， 而本文的滤波算法受到低于地表的噪声点影响 ， 因此噪声点精确剔除是今后的研究重点 。相比 MABEL 试验 ， ICESat-2 的飞行高度更高 、返回能量更弱 ， 因此点云密度很可能比飞行试验更低 ，噪声水平比飞行试验更高 。因而 ， ICESat-2 数据处理对点云去噪算法的要求更高 。同时由于 I CESat-2轨迹长 、覆盖密度高 ， 数据量巨大 ， 所以需要深入研究海量剖面点云的快速滤波方法 ， 为 ICESat-2 在全球变化以及其他领域的应

43、用提供技术支持 。志 谢 感谢美国 NASA 提供 MABEL 试验的数据产品 ; 感谢科罗拉多州立大学激光雷达生态应用中心 ( Center for the Ecological Applications of Li-DA) 提供数据浏览工具 MabelBrowser。参考文献 ( eferences)Abdalati W， Zwally H J， Bindschadler ， Csatho B， Farrell S L， Fric-ker H A， Harding D， Kwok ， Lefsky M and Markus T 2010 TheICESat-2 laser altimetr

44、y mission Proceedings of the IEEE， 98( 5) :735 751 DOI: 101109/JPOC20092034765Anthony W Y， Stephen M A， Li S X， Shaw G B， Seas A， Dowdye E，Troupaki E， Liivap， Poulios D and Mascetti K 2010 Spacelasertransmitter development for ICESat-2 mission SPIE LASE Interna-tional Society for Optics and Photonic

45、s San Francisco: SPIE:757809 757809 11 DOI: 101117/12843342Awadallah M， Ghannam S， Abbott A L and Ghanem A 2013 Activecontour models for extracting ground and forest canopy curves fromdiscrete laser altimeter data/ / Proceedings: 13th International Con-ference on LiDA Applications for Assessing Fore

46、st EcosystemsBeijing: SilviLaser: 129 136Axelsson P 2000 DEM generation from laser scanner data using adap-tive TIN models International Archives of Photogrammetry and e-mote Sensing， 33( B4/1; PAT 4) : 111 118Brunt K M， Neumann T， Markus T， Cook W， Hart W， Webb C， Dimar-zio J， Hancock D， Lee J and

47、Bhardwaj S 2011 MABEL photon-counting altimetry data for ICESat-2 simulations / /AGU Fall Meet-ing， San Francisco: AGU: 1: 0462Brunt K M， Neumann T A， Walsh K M and Markus T 2014 Determi-nation of local slope on the greenland ice sheet using a multi-beamphoton-counting Lidar in preparation for the I

48、CESat-2 mission IEEEGeoscience and emote Sensing Letters， 11( 5) : 935 939 DOI:101109/LGL20132282217Dong L X 2008 Estimation of the forest canopy height and abovegroundbiomass based on multiplicate remote sensing in Three Gorges Bei-jing: Institute of emote Sensing Applications: 18 20 ( 董 立新 2008 基于

49、多源遥感数据的三峡库区森林冠层高度与生物量估算方法研究 北京 : 中国科学院遥感应用研究所 : 18 20)Evans T 2013 Optical development system lab alignment solutions forthe ICESat-2 ATLAS instrument / / Aerospace Conference Big Sky，MT: IEEE: 1 11 DOI: 101109/AEO20136496939Harding D 2009 Pulsed laser altimeter ranging techniques and implica-tions for terrain mapping / / Topographic Laser anging and Scan-ning Principles and Processing UK: Toth CC Press: 173 195 DOI: 1