基于无人机激光雷达遥感的亚热带常绿阔叶林群落垂直结构分析

解宇阳,王 彬,姚 扬,杨 琅,高 媛,3,张志明,*，林露湘

1 云南大学生态学与环境学院暨云南省高原山地生态与退化环境修复重点实验室, 昆明 650091 2 北京大学深圳研究生院城市规划与设计学院, 深圳 518055 3 西南林业大学环境修复与健康研究院, 昆明 650224 4 中国科学院西双版纳热带植物园热带森林生态学重点实验室, 昆明 650223

植物群落结构,包括不同物种的高度和种群密度及其在不同生境中的变化,是植被生态学的基本研究内容[1],也是当前全球气候变化生态研究和大尺度生态系统功能及生物多样性评估的基础信息[2- 4],例如：Schut等[2]在全球生物多样性热点地区探究了通过植被结构变化对气候变化预测的方法；Vogeler等[3]以及Zellweger等[5]实验均表明了植被结构研究在预测森林鸟类物种丰富度中起到重要促进作用；Shugart等[6]利用遥感技术在研究植被结构基础上研究森林碳循环过程。Simonson 等[1]综述了生态学文献中最常用的15种植物群落结构指标包括8种表征植被的垂直结构,在此研究中,水平结构反映景观尺度的斑块特征,如不同类型植被覆盖率、植被类型斑块特征等；垂直结构反映植被竖直空间的变化以及单株树木的数量及形态,如植被冠层高度与树木分布格局[1]。植被冠层高度反映了植物群落在竖直方向占据空间的大小,树木分布格局是生物群落中各种生物与环境因素长期相互作用的结果,直接反映了地表植被的覆盖情况及群落形成与维持机制[7],二者均是影响生物多样性、碳储量及生物量分配等重要因素[8- 9],对于衡量森林生态系统功能方面着重要的意义[10- 11],植被高度在生物量计算的异速增长模型中[12- 13]是重要形态学变量[14- 15],而结合树木分布格局的研究有利于单木尺度平均生物量以及生物量分布格局的测算。

常绿阔叶林在我国是一种占有较大的区域主要森林类型[16],相对于针叶林有乔木层空间体积大、树冠紧凑浓密等结构特征[17]。常绿阔叶林群落结构尤其是垂直结构是进行林下更新等研究的基础[18]。然而当前关于群落结构的研究在常绿阔叶林研究中相对群落演替、更新等方面较少[18],且对于群落结构参数的测量方法以基于传统方法为主。在传统的植被生态研究中,群落结构指标的测定依赖于对植物个体的逐一测量,因此只能设计各种抽样方法来获取样本值,群落冠层高度通常是用测高仪测量几棵或十几棵树的高度求平均值获得[19],树木分布格局的测量基于单棵树木位置的确定,需对样地内树木位置使用GPS一一获取地理坐标[20- 21],因此如上传统方法在获取相对较大尺度或人力难以到达的范围内准确的植物群落结构参数是一项耗时而低效的工作[22],并且传统手工无法进行重复性取样使得测量结果包含空间异质性信息,因此获取的样点或样方等数据结果尺度推演到景观以上水平往往不可靠[23]。随着卫星遥感影像的发展,遥感测量技术为获取大尺度范围的植被结构提供了一个全新的解决途径[24- 25],近些年的研究中,Tao等[26]利用全球卫星数据与实地样方测量,结合气候指标评估了全球植被冠层高度分布模式；Zhang等[9]利用遥感进行了全球最大森林高度分布制图；Crowther等[27]则在全球尺度上进行了树木分布密度的精确制图。而在较小尺度的森林生态系统监测工作中,传统的卫星遥感影像的精度往往难以满足,且基于光学原理的影像数据无法直观反映样地三维特征尤其是林下的地形特征,多种垂直结构需要通过植被指数的计算进行反演,而在此过程中受拍摄时光线条件影响较大,难以反映群落的细节。近年来激光雷达(Light Detection and Ranging,LiDAR)技术快速发展,其作为主动遥感的方法有效穿透森林,相对于传统的光谱遥感能够更好的提取植被三维特征[28],在植被结构参数提取中得到越来越多的应用[29]。其中无人机载激光雷达(UAV-LiDAR)结合了无人机低空摄影测量遥感的技术获取影像分辨率高、重叠率大、相幅小、姿态角大、廉价且不受云雾干扰等特点[30- 33],在获取森林植被冠层高分辨率结构信息方面体现出巨大优势[28]。国内外较小尺度上进行的植被高度与树木密度的研究很多是基于无人机激光雷达开展[34- 36]。但研究研究对象多是针叶林植被[37- 38]。针叶林树冠稀疏,利于激光穿透获取下层植被和地面点云数据,可以更好地还原地形进而提取植被高度等信息,且针叶树冠较阔叶树冠形状规则,更有利于单木提取；而阔叶林情况相反,因此基于激光雷达提取亚热带阔叶林群落结构参数有着更大的挑战性[39]。

本实验利用无人机搭载激光雷达设备对实验样地内的森林群落结构参数包括植被高度以及树木分布格局进行测量提取,对照样地森林群落的地面实测值进行遥感植被参数精度验证,目的在于评价现有的无人机载激光雷达遥感技术及植被结构参数提取方法应用于常绿阔叶林植被结构精细测量和可行性和可靠性,并进行不确定性分析。

1 研究数据与方法

1.1 研究区域概况

本研究选择云南省哀牢山国家级自然保护区范围内的亚热带中山湿性常绿阔叶林作为研究区。云南哀牢山国家级自然保护区纬度范围24°45′35″—23°59′34″N,经度范围100°53′26″—101°29′22″E,总面积约670km2,最高海拔约3100m,一般在2000m以上,高差达2000m以上。年平均降雨量为 1931 mm,雨季降雨量占到年降雨量的85%。年平均蒸发量为1485 mm,年平均温度为11.3℃[40]。有山顶苔藓矮林、中山湿性、半湿润和季风常绿阔叶林4种常绿阔叶林,土壤为山地黄棕壤和红壤[41]。本研究以中科院版纳植物园哀牢山森林生态站为中心选取3个100m×100m的正方形样(图1)地作为研究区域进行无人机及实地数据采集,样地选取依据以下条件随机选取(1)无人机及人力可以到达并进行相应的数据采集；(2)坡度坡向变化复杂以反映真实山区地形特征；(3)植被组成、郁闭度相似；(4)人为干扰较少；(5)样地间距离1000m以上,相互间无影响。样地A地理位置24°32′45.86″—24°32′49.04″N,101°01′11.48″—101°01′15.08″E,海拔2508—2522m；样地B地理位置24°32′33.52″—24°32′36.89″N,101°01′47.12″—101°01′50.71″E,海拔2539—2560m；样地C地理位置24°33′05.42″—24°33′08.71″N,101°01′39.35″—101°01′42.82″E,海拔2586—2598m。群落类型为典型中山湿性常绿阔叶林,优势树种为木果柯(Lithocarpusxylocarpus)、硬斗柯(Lithocarpushancei)、变色锥 (Castanopsisrufescens)等壳斗科(Fagaceae)高大乔木[16],乔木层高度平均20—25m,林相完整。

图1 试验地区与样地位置图Fig.1 The location of the test area and plots

1.2 实验设备概况

本次实验选用数字绿土GV1500 II型号8旋翼无人机作为飞行平台,选取Li-Air无人机激光雷达系统,采用 World Geodetic System 1984(WGS84)地理坐标以及横轴墨卡托投影(UTM),使用华测RTK系统定位,激光扫描头不断地记录不同方向的测距数据,从而得出以扫描中心为基准的三维坐标信息,发射激光波长为905 nm,测距范围为100m,但在60m范围内最佳。

1.3 数据采集

1.3.1无人机激光雷达数据采集

完成一次无人机近地面低空摄影作业,其步骤是：无人机航线设计规划——航线上传给无人机并验证——无人机搭载传感器自动按航线执行飞行任务——返航降落并收集数据。Li-Air无人机激光雷达系统激光扫描头测距在60m范围内时获取的激光点云质量最佳。航线设定的高度是相对于地面固定基站的相对高度,在地形起伏较大的山区,为保证激光扫描头到地物的垂直距离始终≤60m,需设置变高飞行航线。通过无人机地面站软件基于Google Map框出1hm2样地范围,自动生成航线水平轨迹与航线拐点处若干目标点,目标点可导出为带有经纬度信息的shp点文件。另外需获取样地遥感图像,本实验利用无人机搭载MicaSenseRedEdge多光谱相机,设置全程飞行高度200m,旁向重叠度65%,获取样地多光谱影像(图1),利用AgisoftPhotoscan软件立体像对生成正射影像并采样内插生成高精度的数字表面模型(Digital Surface Model, DSM)。将目标点shp文件与栅格数据DSM导入ArcGIS,使用Extract Values to Points算法按目标点位置提取对应DSM像元值即海拔数值,该数值减去地面基站海拔高度值加上60m即为该目标点需设置的高度值。在地面站软件上进行设置即可生成变高飞行航线。无人机搭载Li-Air激光雷达系统执行该航线,获取激光雷达原始数据。

1.3.2样地实地数据采集

使用与激光雷达相同一套华测RTK移动站记录3个样地内所有上层乔木的位置信息,使之与激光雷达数据统一在相同的地理及投影坐标系中。使用用激光测高仪在每个样地中实地随机测量40株上层乔木的高度进行记录。

1.4 激光雷达数据获取与去噪

分别从激光雷达存储控制单元与地面基站下载当次作业原始记录数据以及基站移动站数据(POS数据)。耦合结算完毕即获得样地的激光雷达点云数据(las.格式)。此时的激光雷达点云数据由于受到空中云、鸟等干扰而产生中心点明显高于周围点的平均高程或因为多路径反射而产生的明显比实际点位低的噪点。对数据中每一个激光点搜索指定邻域点个数的邻点,计算该点到邻点的距离平均值并计算所有平均值距离的中值和标准差,如果某点的平均值距离大于最大距离(中值+标准差倍数*标准差),则将该点分类至噪点类别。分类后仅导出未分类类别,即可获得样地无噪点点云数据(图2)。使用Terrasolid软件拓展模块TerraScan的classify air/low points算法可实现上述噪点去除。

图2 样地激光雷达点云数据Fig.2 LiDAR point cloud data of each plot

1.5 植被树冠高度提取与精度检验

1.5.1树冠高度模型提取

提取植被冠层高度的关键是冠层高度模型(canopy height model,CHM)[4],由于样地地表覆盖地物可以确定只有植被,因此可以直接将由冠层上表面的点生成数字表面模型(DSM)与利用地面点生成的数字地形模型(Digital Terrain Model, DTM)做差得到,即DSM-DTM=CHM。DSM可直接基于去噪后的未分类点云导出最大值矢量格式灰度图获得,即得到的栅格数据每个像元值为该位置上最大高程激光点海拔值。DTM获取的关键在于地面点的确定,根据激光雷达激光反射的特点确定点云最低点构成数据全部地面点。寻找地面点的算法为反复建立地表不规则三角网模型[42],首先选取少量点,建立初始稀疏三角网模型,初始模型的三角形大多数低于实际地面,只有较高点接触到地表。算法通过设定的反复参数(Iteration angel和Iteration distance)反复加入新的激光点(空间插值)开始向上扩建模型,每个加入的点使模型更加贴近地表。反复参数决定一个点有多近才能被纳入三角平面[43]。本实验设置Iteration angel即选中点和三角网格的最近顶点的连线与该三角网格平面的最大夹角为10°,Iteration distance即点与三角形的最大垂直距离为1.5m,地表所允许的大坡度90°,基于上述参数将地表点提取出并生成阴影面,导出矢量格式灰度图,即得到样地DTM栅格数据。将DSM与DTM 数据相减得到这两个栅格数据的差值即CHM数据,每个像元值即该处植被冠层高度。上述DSM以及DTM的提取使用Terrasolid软件拓展模块TerraScan；二者做差得到CHM在ArcGIS中使用栅格计算器工具实现。

1.5.2精度检验

对提取到的植被高度进行精度检验。将每个样地实测的40棵树木信息(包括地理位置信息与高度信息)转成点矢量格式导入ArcGIS,赋予每个实测位置点对应CHM栅格值。精度计算公式为公式：

(1)

式中实测值为用激光测高仪实地测的树木高度,估测值为实测树木对应位置上CHM值。计算均方根误差(RMSE)衡量观测值同真值之间偏差,并对对实测值与估测值进行线性回归分析验证相关性。

1.6 树木分布格局计算与精度检验

1.6.1树木位置提取及精度检验

基于单木树冠分割算法提取树木位置。本实验用局部最大值法搜索树冠顶点[44- 45],对于阔叶树来说,其树冠近似呈现中心高四周低的圆拱,因此搜索到的树冠中心最高点可以用来定义树木位置[46]。对之前提取到的CHM数据过滤从而提取CHM值最大的点视为树冠的位置,提取间隔可通过设置最小冠幅x来确定,即设置最大值提取间隔,比较每一个像元与以其为中心周围x距离范围内的像元,提取出值最大的像元即为树冠顶点。本实验根据实地观测设置最小冠幅即树冠半径为1m,在LiDAR 360(试用版)软件中实现上述树木位置提取的操作。

对树木位置提取结果进行精度检验。本实验的树木位置提取精度检验采取样地尺度与单木尺度的方法[47- 48],样地尺度检验从样地整体角度出发,不考虑提取单木位置和实测单木位置的对应关系,所需指标为单木探测百分比(Detection Percentage,DP),即样地内提取出的单木数量占实测单木总数量的比例,公式表达为(2)：

(2)

式中,Nr为实测样地内单木总数量；Nd为LiDAR探测到的单木总数量。

单木尺度位置检验所需指标包括1∶1对应关系单木数量(N1∶1,本实验中,实测得到的单木与其在其实测位置的2m缓冲区内、所探测到的激光雷达提取出的直线距离最近的单木位置之间形成的关系被称为“1∶1对应关系”,这些实测出的以及激光雷达提取出的单木被称为“1∶1对应关系单木”)、用户精度(User′s accuracy,UA)以及生产者精度(Producer′s accuracy,PA)[48]。其中后两者的计算公式如公式(3)(4)：

(3)

(4)

计算得到距离实测点2m范围内的激光雷达估测点；再基于后者建立半径2m的buffer裁剪前者,剔除未在2m范围内探测到估测点的实测点；比较现阶段提出的两种点数量,基于数量较少者,计算得到其每个点到较多者的最近点以及二者的直线距离,找到多对1的情况,即多个数量较少点距离最近的为同一个数量较多点,挑出较远者分配与其第二接近且在2m范围内的较多点对应,如2m内无其他点则剔除。最终获得的数量较少者数量即为1∶1对应关系单木数量。通过上述方法得到提取结果的单木探测百分比、用户精度以及生产者精度等信息,并计算1∶1单木间平均距离。上述精度检验过程在ArcGIS中进行。

1.6.2树木分布格局计算

采用Clark-Evans最近邻体指数法计算空间格局,即通过相邻最近树木距离平均值与随机分布下所期望的平均距离之比确定树木分布的聚集程度[49]。同时为了减少斑块的大小和形状对空间格局分析的影响,考虑到最近邻体法在GIS中运用的有效性与实用性,故选用Fǜldner提出的最近邻体法修正公式计算[50- 51],同时考虑到样地或斑块的形状与边缘效应[52],提高分析精度,公式如(5)所示：

(5)

式中,CE表示Clark-Evans指数,rA为样地或斑块中每株单木个体与其最近邻体间距离的平均值(m),rE为该样地或斑块内所有个体随机分布时rA的预期值(m),ri表示第i株树木个体与其最近邻个体间距离(m)；N为该样地或斑块内单木个体总数,A为该样地或斑块面积,P为该样地或斑块周长。rA与rE的偏离程度可通过正态分布检验如公式(6)：

(6)

(7)

式中,σ为标准差,ρ为公式(7)即斑块内树木个体密度,其中N为样地或斑块内树木个体数量,A为该样地或斑块面积。按照正态分布检验原则实际如果μ<1.96(即显著水平为0. 05时的临界值),从统计意义上可认为CE=1, 单木分布格局判断为随机分布。如果μ>1.96,CE<1可认为实测CE值显著小于1, 判断为聚集分布；CE=1,判断为随机分布；CE>1,判断为均匀分布。对激光雷达提取结果与实测结果分别进行上述操作,后者作为前者对照。上述分析过程在ArcGIS中进行。

2 结果分析

2.1 植被高度提取结果与精度检验

基于上述方法提取到的样地CHM如图3所示,树高实测值与对应CHM值线性回归分析如图4所示,精度分析如表1所示。结果显示平均精度达到94%以上,均方根误差总体1.28m,3块样地线型回归系数平均0.92以上,P均<0.01,均有显著统计学意义,R2分别为0.905、0.909、0.883,总体拟合度较高,激光雷达提取的植被高度值与实际值接近,在实际研究中可代替实地测量应用。

图3 冠层高度模型(CHM)及树木位置点Fig.3 Canopy Height Model (CHM) and Individual trees location

2.2 树木位置提取结果精度检验

LiDAR数据提取到的树木位置点如图3所示,树木位置提取精度评价结果见表2。在样地尺度上,各样地探测百分比都在80%以上,平均为86.01%,样地C最高,接近90%,总体精度较高,但漏测现象普遍存在,即LiDAR探测到的树木普遍少于实际的树在单木尺度上,由于同一样地内的LiDAR探测单木树木数量均小于实测单木总数量,因此生产者精度均低于用户精度,前者大致在60%—70%之间,后者大致在70%—80%之间,漏测导致了一定的误差。对于提取到的1∶1对应关系单木,平均偏差距离在1m左右,标准差约0.4m,远小于单木平均冠幅半径,因此多数1∶1对应关系单木可看做实际同一单木。

表1 各样地激光雷达植被高度提取结果精度分析

CHM： Canopy height model; LiDAR: Light detection and ranging;

表2 树木位置探测精度评价

Nr：实测样地内单木总数量Total number of individual trees in the plot；Nd：LiDAR探测到的单木总数量Total number of individual trees in the plot；DP：单木探测百分比Detection Percentage；N1∶1：1∶1对应关系单木数量Number of 1∶1 correspondent trees；PA：生产者精度Producer′s accuracy;UA:用户精度User′s accuracy；d1∶1：每一对1∶1对应关系单木平均距离Mean of distance of each pair of 1∶1 correspondent trees

2.3 树木分布格局计算结果

基于修正Clark-Evans最近邻体指数法计算得各样地分布格局及相关参数如表3所示。由于LiDAR的漏测情况导致树木密度测量值略小于实际值。样地实测显示3个样地内上层乔木均呈聚集分布,而基于LiDAR的测量结果为随机或均匀分布。

表3 各样地树木分布格局

AD：聚集分布Aggregation distribution；RD：随机分布random distribution；UD：均匀分布 uniform distribution

3 讨论

目前利用遥感技术研究森林植被指标的方法有很多种,然而传统的光学遥感技术获取森林三维结构的能力有限[53]。激光雷达作为一种主动遥感技术可以有效穿透森林,结合无人机系统发展对空间生态学研究的推动[54],在获取森林垂直结构参数方面有着无法比拟的优势[28]。也有基于无人机获取多光谱数据或可见光数据进行植被高度和单木提取的案例,王彬等[19]通过无人机搭载数码相机获取的人工针叶林数码相片即可见光数据,结合地面设定的坐标控制点赋予影响地理坐标及投影坐标,参考R、G、B、亮度、色调等像元信息分离地面像元与植被像元,并选取提出地面像元点内插生成地表面高程数据以获得地形信息,进一步处理获取树高及位置,这种方法成本较低且操作方便,适合于针叶林或稀树草原等生境,但在植被覆盖率很大、郁闭度高的阔叶林里,难以通过数码图像或多光谱影像直接分辨提取地面点,且受树阴以及山体阴影等影像,更难以直接通过光谱信息进行非植被点提取。由此便体现出激光雷达点云数据的优势,其原理通过发射激光接触到表面后折返来计算目标与自己的相对距离[55],每一个激光点都反应了地物相对真实的位置信息,且系统自带内置的卫星定位系统,有相对准确的测定内方位元素,使得获取的数据自带地理坐标及投影坐标信息,相对于传统的无人机遥感影像无需进行几何校正,节省了在样地中设置地面控制点环节,为研究高郁闭度、地形复杂及人力难以到达的样地提供了便捷。

本实验进一步证明了无人机LiDAR在提取植被高度方面的较高的可行性与准确性。树木位置识别是计算树木密度以及更多单木尺度植被参数的基础,与植被高度结合,是计算森林生物量与平均生物量的基础指标。基于CHM局部最大值法探测到的单木数量以及树木密度均略小于实际值,这与样地较高郁闭度有关,基于激光雷达点云数据可进行植被覆盖度的提取[56],用非地面点数量除以该样地激光点总数即可计算得3个样地植被覆盖度分别为：0.807、0.825、0.702,均属中高度郁闭,根据实地目视测算,高郁闭度条件下相邻树间树冠重叠度较大,相邻树木最高点难以区分,给单木分割增加了一定的难度,这是出现漏测现象的主要原因,另外无人机载激光雷达设备在能够获取较大范围内的森林垂直结构尤其是上层树冠信息的同时,由于森林冠层的遮挡阻碍激光穿透形成多次回波,难以精确获取下层植被信息[28],因此部分没有完全暴露的树冠位置难以提取测算。

在自然条件下,由于生境异质性等原因导致乔木层植物呈聚集分布[21],同时徐远杰等[57]实验证明在哀牢山地区自然条件下的中山湿性常绿阔叶林中上层优势树种多为聚集分布,本实验的实地测量对上述研究成果进行了进一步证实,而基于LiDAR进行的树木空间分布格局的判定为随机或均匀分布,较实际情况有所不同,从数据特征角度分析,LiDAR提取的树木位置数据由于是基于局部最大值提取的,实际得到的是树冠顶端,与实测的树木位置相比,基于前文提取的1∶1单木数据来看,有着1m左右的偏差,而在实际的森林里,由于竞争和密度制约会使聚集度下降以致均匀或随机分布[58],推测树冠由于竞争阳光等原因会呈均匀分布格局,因此基于LiDAR提取的数据聚集程度小,偏向均匀或随机分布。另外本次实验将全部乔木作为研究对象,没有进行物种划分,而前人多数的植被空间格局的研究多是针对单个物种,由于种间竞争、资源配置等原因单个物种在空间内呈聚集分布。本实验中不同乔木树种在顶极群落内达到一种稳定的状态,因此形成均匀分布的格局也有其合理性。

本实验尚且存在很多不足需加以改进,如本次获取的激光雷达点云转化为栅格后最高的精度仅为0.5m,相对于无人机近地面遥感所追求的cm级数据有着一定的差距,在未来的研究中,可以通过进一步精细规划航线和采用更高级无人机激光雷达设备等方式提高分辨率。另外精度分析可知基于无人机LiDAR通过局部最大值法提取常绿阔叶林树木位置结果尚且不能真正替代实际测量的数据进行使用,研究应进一步探究不同的适用于高郁闭度的阔叶林单木提取方法,如结合地基激光雷达获取林下植被,基于树干自下而上获取位置信息[8],并结合现有的无人机LiDAR进行多源LiDAR数据融合提取完整的树木。

4 结论

本研究进一步表明无人机载激光雷达遥感技术应用于植被结构精细测量的准确性。选取哀牢山3块1hm2中山湿性常绿阔叶林样地验证激光雷达提取冠层高度及树木位置精度,得到的植被高度数据与实地测量数据对照,二者有极显著相关性,三块样地平均精度在95%左右,实验证明基于激光雷达可在实际研究中代替实际测量提取植被高度数据。基于CHM使用局部最大值法获取树木位置信息,单木探测百分比平均在86%左右,用户精度以及生产者精度平均分别为75.69%和65.15%,主要是由于群落郁闭度及树冠的冠幅和连接度等因素导致部分漏测的现象,1∶1单木间平均距离为1.23m,总体来说树木位置精度较高,但基于LiDAR的树木分布格局显示高大乔木呈均匀或随机分布,与实际情况下的聚集分布不符。在未来的研究中,应进一步研究更为精确的单木提取以及植被高度提取方法,为通过无人机激光雷达测算森林生物量提供更加精准的指标；同时开发新的算法实现基于激光点云数据的植物物种分类,以探究不同种的垂直结构性质。

致谢：感谢中国科学院西双版纳热带植物园哀牢山生态站提供的实验条件,感谢北京数字绿土公司提供的无人机激光雷达系统及数据解译与处理软件。