机载激光雷达树高产品的释义与精度水平评估

付安民,高显连,吴发云,高金萍

(国家林业和草原局调查规划设计院,北京 100714)

树高是最为重要的森林定量遥感观测参数之一,高空间分辨率、高精度的树高遥感产品是进一步开展准确森林资源定量遥感调查和生态环境定量遥感评估的必要数据基础[1]。目前，大范围、高质量的树高遥感产品主要由航空激光雷达观测数据处理得到[2],并且随着近几十年来激光雷达观测和计算机技术的快速发展,航空激光雷达数据的采集和处理能力不断得到提升[3]。激光点云数据已渐渐成为我国森林领域的重要数据源,服务于国家和区域性的森林资源调查、森林碳储量估算、森林生态多样性评估等工作[4-7]。

从激光点云到树高的过程通常需要分别获取林下地形(Digital Terrain Model,DTM)和林冠表面高度(Canopy Surface Model，CSM)产品,两者的差值作为森林冠层高度产品(Canopy Height Model,CHM),因而测量树高的能力受冠层表层和林下地形高度双重精度限制,获取高质量的CHM是一项具有一定挑战性的工作[8-9]。从观测对象上看,森林观测场景复杂多样,不仅树种组成多样,而且林木高低起伏、相互遮挡、枝叶纵横交错,并且耦合了地形要素。“完美”的数据采集不仅要高质量捕捉到林下地形的起伏,同时还要准确刻画林冠大小和高低,才能保障高分辨率林下地形和林冠表层高度产品的准确性和精度。从激光雷达载荷角度上看,不同载荷或同一载荷不同参数设置,观测的准确性和精度指标均会受到影响,会引起观测森林树高的不同。由于森林是非刚体的“软”目标,激光脉冲波束对树冠具有一定的穿透性,直至照射至地面；激光探测器记录透射激光脉冲的回波波形,依据波形振幅分解的点云结果进行地面点分类滤波处理,同步获取林木冠层和林下地形高度。对两者的测量能力取决于激光雷达载荷波形探测器的敏感性及其波形分解算法能力、激光脉冲能量、观测几何条件,以及激光照射地面光斑大小等综合因素[10]。

林冠高度产品CHM并非真正意义上的地面调查树高,其物理意义需要与地面实际树高进行标校释义[8]。之后,通过进一步的地面树高归一化处理,以达到联系林木蓄积、生物量模型,实现森林资源调查的应用目的。同时,受到山区林地复杂的测量条件,激光雷达观测条件差异和处理技术方式等综合因素的影响,CHM计量的准确性及其精度指标在测区空间分布上具有不确定性,有必要建立地面林木参照对象,对产品质量进行系统的检验分析,以促进激光点云采集数据指标和数据处理方法的优化改进,为获取高空间分辨率、高精度的树高遥感产品提供保障。

本文以航空激光雷达森林资源调查数据工程质量控制为目的,设计了航空激光雷达CHM高度产品评估规程方法,旨在标称CHM的“真实”物理意义,分析CHM高度产品的准确性和精度指标水平,检验样地数据与航空激光点云数据的一致性。

1 研究区和数据

1.1 研究区

2018年9月1日—9月27日,在吉林省东部东北虎豹国家公园部分地区(宁安、汪清和穆棱境内),地理坐标(43°19′～44°10′N,129°4′～130°21′E),东西长108km,南北宽105km范围内,开展了总面积 5 600km2的航空飞行遥感森林调查综合实验,如图1所示。实验以复杂山区地形条件下航空森林资源调查为主要目的,同步开展地面森林资源数据调查。

研究区属于温带大陆性季风气候,气候特征表现为春季多风少雨干旱,夏季炎热短促,秋季冷凉降温迅速,冬季寒冷漫长。年降水量变化在450～750mm,降水集中在5—9月。研究区以森林生态系统为主体,主要是温带针阔叶混交林,其中,次生林分布广泛,原生性红松阔叶混交林仅呈零星分布[11]。

注:方框为312幅激光雷达文件数据,红点表示地面森林调查样地。

1.2 数据

1.2.1地面林木调查数据

地面森林调查结合历史森林资源规划调查数据,进行了10个类型:7个纯林(云杉(PiceaasperataMast.)、冷杉(Abiesfabri(Mast.)Craib)、落叶松(Larixgmelinii(Rupr.)Kuzen.)、桦木(Betula)、杨树(PopulusL.)、椴树(TiliatuanSzyszyl.)、栎类(QuercusL.))和3个混交林(针叶混、针阔混、阔叶混)的抽样设计,共计完成439个有效样地调查工作,分布情况如图1所示。样地为圆形,直径大小为30m。样地内胸径大于5cm的所有样木均划为检尺范围,所有检尺样木记录树种、胸径、树高,部分样木记录冠幅、枝下高参数。采用微波、激光测距仪进行树高和冠幅的量测记录,采用胸径尺进行胸径(树干高度1.3m处)测量；相关林木指标如表1所示。

样地中心和单样木水平位置,首先协同卫星定位服务连续运行参照站(Continuously Operating Reference Station,CORS)利用差分技术实现双点绝对定位；然后再架设全站仪利用后方交会测量方法实现单木精密水平定位,名义上样地中心和单样木水平相对精度可优于0.5m,绝对精度优于2m；相关林木定位指标如表1所示。

表1 样地样木野外测量指标

1.2.2航空激光雷达数据

航空激光雷达数据采集利用机载激光雷达扫描系统联合IMU/DGPS辅助航空测量技术开展,具体选用塞斯纳208B飞机搭载RIEGL-VQ-1560i激光雷达载荷,集成惯性导航(Inertial Measurement Unit,IMU)和全球定位导航系统(Global Navigation Satellite System,GNSS),并联合CORS站点同步观测,具体设备及其飞行指标如表2所示。

实际飞行数据达到激光脉冲密度优于10pulse/m2。飞行任务规划为南、北两个测区,飞行高度分别为2 100m和2 300m,飞行速度均为260km/h,飞行航带重叠度优于25%,激光测距仪载荷的参数设置保持一致。南、北两测区分别设置了构架航带,构架线航测区域布设校准面,对飞行航线进行整体校准。

表2 航空激光雷达数据采集设备及其指标

2 原理和技术方法

2.1 样地样木数据处理

样地样木数据处理过程主要分为4个步骤:1)GNSS数据进行差分定位处理；2)协同全站仪后方交会测量数据和GNSS差分结果的单木定位和样地中心定位；3)样地林木统计；4)建立森林资源样地和样木数据库。其中GNSS绝对定位精度优于2m；全站仪单木相对定位精度优于0.5m。平面坐标为CGCS 2000坐标系,高程采用1985国家高程基准,投影采用高斯克里格,3度分带,中央子午线129°。

其中样木的水平定位需要特别注意,开展林区定位测量,GNSS接收机接收卫星信号易受到地形、林木遮挡,经常发生卫星信号多路径和无法直接接收情况,定位精度受到影响,如图2所示。通常需要延长GNSS接收机的观测时间达到15min以上,才能达到绝对定位精度优于2m的需求。此外,事后GNSS差分软件报告精度与实际位置精度也没有相关性,无法通过软件报告的精度水平,推断实际的定位水平[12]。因而激光点云冠层高度模型数据CHM与地面样木进行高度对比时,需要进行水平一致性检验,并对偏移明显的样地进行修正,使得两者的位置匹配精度优于2m。

图2 GNSS卫星接收信号多路径效应和无法直接

2.2 激光雷达数据及DTM CHM产品

激光雷达数据处理主要分为4个步骤:1)利用POSPac软件将GNSS和IMU惯导数据,联合地面参照站数据进行差分定位处理,逐航带进行定位定向精结算；2)利用RiPROCESS软件将精解算姿态和定位数据赋给每个激光脉冲数据,并根据架构航线对各航带激光雷达数据进行航带间拼接和整体平差；3)剔除明显的噪声点；4)区域分幅,对激光点云数据进行滤波和分类,区分地面点和非地面点。东北虎豹国家公园机载激光雷达数据2,3级产品融合一体,共计312幅文件,总数据量2.7TB,如图1和表3所示。

在激光雷达2,3级数据的基础上,首先,利用地面分类点,构建TIN网空间插值生成1m空间分辨率的DTM数据；其次,利用首次回波点云(First Return)数据,根据1m网格高程最大值统计方法,以非空间插值方式,获取CSM数据；最后,两者差值获取1m空间分辨率的CHM数据。

2.3 CHM质量分析方法

CHM质量分析方法主要有6个步骤。1)确定样木空间尺度大小,假定单木冠幅直径不小于5m。在此假设条件下,以样木水平定位坐标为中心,建立2m半径的圆形缓冲区,用来分析样木彼此的空间临近关系,执行空间叠加运算分析。2)优化选择地面参照对比样木:为避免林木高低起伏,相互遮挡影响,运用空间拓扑分析和局部树高最高原则,逐个样地选择具备较好观测天窗的上层林木,且确保周边没有其他林木遮挡。3)空间叠加分析:将地面样木数据与CHM产品空间叠加,逐个样木获取缓冲区内高度最大值CHMmax,赋予每株样木；建立CHM树高、地面调查树高的双样本集。4)从样木总体看CHM树高:根据双样本散点图和差值直方图统计量,反复迭代去除粗差(均值μ及2倍标准差2σ为参变量)；利用优化样本总体确立通用型CHM树高归一化模型,对标CHM高度“真实”物理意义,评估其准确性。5)从分树种组样木总体看CHM树高:根据双样本散点图和差值直方图统计量,查看各树种组CHM高度“真实”物理意义和准确性分异情况。6)评估CHM产品质量:制定精度水平分级指标,统计抽样参照对比样木的CHM测量精度分级情况,指示CHM质量状况。

表3 激光点云数据产品分级表

图3 CHM产品生产与精度评价流程图

3 结果与分析

3.1 CHM树高准确性和精度指标分析

从地面林木调查数据中选取上层林优势木作为“真实”数据,用于CHM树高准确性和精度的对标分析。首先,从地面每株林木定位精度优于2m条件出发,构建每一株样木的2m半径圆形缓冲区；其次,利用空间分析工具确定局部优势的高大乔木,即挑选的每株林木在与其缓冲区交叠的所有邻近林木中具有优势高度；剔除周边低矮林木。图4显示了初步优化后在中密度林和高密度林中的典型抽样结果。最终,从439个样地的全体25 467株林木中,挑选了参照样木7 873株。

将确立的参照样木按照缓冲区分析范围提取最高值CHMmax作为机载激光雷达树高,建立双样本对比数据,图5(a)和(b)显示了两者的散点图和差值直方图(参照样木树高减去CHMmax)。从图5(a)散点图上看,两者存在明显的线性关系,R2=0.68,同时中低林木观测存在显著的CHM树高大于地面参照样木的现象,大多是被临近高大林木遮挡导致的,如图4所示。从图5(b)树高差值直方图上看,在-2.5～4m的区间范围呈现明显的高斯峰分布；同时还存有700余株的CHM树高显著大于地面参照样木现象(5m以上)。综上,CHM树高和地面参照样木双样本呈现显著的线性相关性,且CHM树高略低于地面参照树高,均值约为1m左右；但是也存在一定数量的粗差样本对。

注:黑色数字表示地面调查树高,红色数字表示地面树高与CHMmax的差值,高亮圆圈表示受邻近高大林木影响导致的粗差样本。

图5 地面参照样木树高与相应2m半径圆形缓冲区CHMmax的散点图和差值直方图(迭代去粗差前后对比)

去除粗差样本,优化上层林木确立“真实”数据样本集,用于树高物理意义的对标分析。以μ±2σ为限值(样本均值μ和标准差σ),经4次去粗差迭代,去除粗差样本。重新获取的优化上层林木样本总计5 856株,如,图5(c)和图5(d)显示了总体散点图和树高差值统计直方图。两者线性相关系数R2=0.97,树高差均值μ为0.7m,标准差σ为±1.0m。表4显示了各个树种(组)的专题统计结果(1)吉林省林业厅.吉林省立木材积、出材率表.2016.。从各个树种(组)的树高差均值μ看,均呈现CHM树高比上层林木树高要低的现象,偏差范围0.3～1.3m不等,其中偏差较大的为人工针叶林Ⅰ的1.0m、人工针叶林Ⅲ的1.3m；从各树种(组)树高差变量标准差上看,大小范围±0.7～±1.2m不等。以上分析表明,CHM树高与地面真实树高存在系统偏差,各个树种略有差异。

表4 CHM树高和精度标校分析表

3.2 CHM树高产品精度评价

CHM树高产品精度评价就是筛选适合的上层林木作为“真实”树高,开展对比分析,定量化反映产品的精度水平及其空间分异情况。其中,地面调查林木的单株水平定位精度优于2m；同时,地面树高测量精度受林地复杂观测条件限制,通常误差范围允许最大达到“真实”树高的10%。这样的条件表明,CHM数据总体精度指标不能由个别树高决定,而是应从样地中抽取一定数量上层林参照样木,通过数据统计,更为可靠地指示CHM数据准确性和精度误差；在局部样地尺度,也应形成一定数量团组抽样模式,达到指示样地单元精度水平的目的。

从图4可以看出,2m半径的圆形缓冲区内有时会涵盖邻近高大林木的树冠,因而CHMmax的算法设计导致了有些林木树高提取不正确。显然,需进一步优化上层林样木,主要采取以下两个步骤:一是去除样地外围高大林木的遮挡影响；二是提高上层林木缓冲区窗口尺寸,尽量剔除周边低矮林木。具体设置参数为:依次采用2,3,4,5m半径的圆形缓冲区窗口；设置原样地范围边界向内缩减5m。随着缓冲区大小/样地半径变化,从表5看出,地面参照林木的株数依次为7 873株(5),3 366株(4),1 601株(3),1 017株(2),714株(1),逐步快速递减。这与保障每个样地保持一定数量抽检样木如5株以上的目标矛盾。为平衡样地尺度数量问题,在方案(5)基础上,每个样地抽取5株最高林木,作为补充抽样方案(6)。

3.1章节中初步明确了正常情况下大部分地面参照林木与CHM树高的总体精度水平,即CHM树高与地面参照林木树高存在系统误差μ=0.7m,标准方差σ=±1m。据此采用标准差分级方法,以算术平均值μ作为中间0点,以±σ,±2σ,±3σ,±3σ以上/下为分界点,进行分组计数统计,从而指示CHM的总体精度水平和显著异常数据。从表5看,方案(6)最为有效的抑制了林木遮挡影响,-3σ以下的分组占比最低,仅为0.7%,可作为评价CHM树高产品精度水平的最优方案。由此得到示范区CHM产品的精度水平状况:其中0～±σ的百分比为57.6%,0～±2σ的百分比为79.4%,-2σ以下占比2.0%,+2σ以上占比18.6%。图6具体展示了表5中抽样方案(5)和(6)的前后样地尺度抽检样木的空间分布和高度误差的对比情况。方案(6)避免了过多粗差样木的干扰,集中通过5株上层林木准确反映了样地尺度的树高误差情况。

表5 地面参照林木抽样方案及其标准差分级统计表

(a)样地参照样木全部合格情况

(b)样地参照样木全部不合格情况

(c)样地参照样木仅部分合格情况

注:黑色数字表示地面调查树高,红色数字表示地面树高与CHMmax的差值。假定CHM树高与真实树高差值介于μ～±2σ为合格抽检样木,

依据表5抽样方案(6)精度指标分级统计结果,+2σ以上占比为18.6%。说明CHM产品主要问题是与实测树高比较CHM树高偏低。虽然激光雷达测量树高的理论十分简单,但是树高的准确性和精度受到仪器性能、采集规范和数据处理诸多因素的影响,十分复杂,具有一定的不确定性。CHM树高误差Etree的3个方面来源[8]分别是:1)冠顶高度偏低误差Etop,其程度大下是林木自身形态如树冠大小和形状、枝叶疏密与激光脉冲能量、脉宽、传感器灵敏度等相互共同作用的结果；2)林下地形DTM数据高估误差EDTM,其程度与采集过程中的地面点质量和地面点滤波处理的有效性密切相关；3)林木自身倾斜导致的误差ELean,其程度与林木倾斜角度和局地地形坡度大小相关。除了以上误差因素外,地面调查数据的错误也参杂其中,如,树高记录错误、林木定位错误等。上述误差和错误可依据指标分级的指示,逐个样地进行人工目视判读,识别误差来源和错误情况如图6所示。

4 结论

航空激光雷达数据及其产品已成为森林管理和森林动态监测的重要技术手段,但是其精度水平目前主要以平坦区域或明显的人工地物控制点进行评价；山区林地由于复杂地形和森林测量条件成为精度评估的“盲区”。随着航空激光雷达在森林调查和管理工作中的广泛应用,大面积林区CHM数据精度指标缺失的问题日渐突出。

本研究从现有地面林木水平定位精度优于2m、测高精度为林木“真实”树高10%的实际工作条件出发,利用空间分析方法和统计分析策略,从工程应用角度制定了CHM树高产品的检验规程。在我国东北虎豹国家公园航空激光雷达森林资源调查示范区的应用结果表明,该方法准确、可靠地反映了地面调查样地、样木与CHM高度数据的空间一致性情况:一是验证了山地林区CHM产品精度水平,指示出偏差异常情况,从而为查找数据采集和处理过程中存在问题提供线索；二是能够检验地面调查样地的样木定位和高度测量的准确性。