基于激光雷达的大别山毛竹冠层高度反演研究

李玲娜，孙 悦，万芳芳，徐小牛，黄兴召

(安徽农业大学林学与园林学院，安徽 合肥 230031)

毛竹(Phyllostachysedulis)适应性强、生长快、经济价值高，也是我国竹类资源中分布最广的竹种。现有的毛竹林主要属于笋材兼用林，是山村林农重要的经济支柱[1]。根据第九次森林资源清查数据，我国竹林面积为6.41×106hm2，其中毛竹林面积4.67×106hm2，占竹林面积的72.9%；毛竹株数1.41×1011株，胸径主要分布在7～11 cm，占61.9%[2-4]。如何进一步精准获取毛竹林的冠层高度，是掌握毛竹资源结构特征及可持续管理利用的关键。

激光雷达是近年来发展较快的一种主动遥感技术，为准确提取森林结构参数提供了有力的技术支撑[5-6]。搭载在有人飞机或无人机上的小光斑激光雷达统称为机载激光雷达(Airborne LiDAR)[7]。利用机载激光雷达数据生成数字高程模型(Digital Elevation Model，DEM)和数字表面模型(Digital Surface Model，DSM)，计算出冠层高度模型(Canopy Height Model，CHM)[8]。机载激光雷达获取高精度的点云数据，能够更加精确的估算毛竹冠层高度[9]，但是机载激光雷达做不到大尺度连续的观测，需要结合其他大尺度连续的数据共同反演毛竹冠层高度，且机载激光雷达数据获取成本大，很难进行大规模的推广[10]。星载激光雷达以卫星为搭载平台，运行轨道高，观测范围广，全天时全天候对地进行观测[11]。高精度的机载激光雷达数据和大尺度的星载激光雷达数据相结合，进一步提升了该方法的应用范围[12]。

全球生态系统动态监测(Global Ecosystem Dynamics Investigation，GEDI)是目前唯一专门为测量植被的垂直结构而设计的全波形激光雷达，其核心目标是估算热带和温带森林的地上生物量，模拟地表碳收支，了解植被结构对生物多样性的影响[13-15]。近年来，有关GEDI反演森林冠层高度已有较多研究，但冠层高度受树种、气候环境、地形等因子的影响，导致模型预估效果存在差异[16-18]。因此，利用GEDI数据，开展不同树种的冠层高度反演研究，是提升模型预估精度的基础[19-22]。为此，以毛竹为研究对象，建立GEDI反演模型，探讨不同因子对模型的影响，并找出影响模型预估的关键因子，以期为毛竹林的结构监测和可持续利用提供技术支撑。

1 试验地概况与研究方法

1.1 试验地概括

实验选择安徽省金寨县关庙乡竹园村(115°25′E，31°30′N)为研究地，该地平均海拔450 m，属于北亚热带季风气候，年平均温度14.6 ℃，极端高温和低温分别为35 ℃和-14 ℃。地带性植被为常绿阔叶林，毛竹林为该地主要的森林植被类型之一。土壤主要是黄棕壤，土层厚度在0.5～1.2 m，土壤pH值在5.5～6.5之间。2021年7月，选择代表性立地条件及竹林结构，建立25 m×25 m的标准样地28块(见图1)。在样地内进行每木检尺，利用超声波树高测定仪(LD6172)测定树高，作为模型预估效果验证数据。

图1 毛竹林临时样地分布图Fig.1 The distribution of Moso bamboo sampling plots

1.2 研究数据获取

1.2.1星载激光雷达数据 GEDI星载全波形激光雷达数据于2021年7月在NASA/USGS (National Aeronautics and Space Administration/United States Geological Survey)土地处理分布式活动档案中心发布，GEDI二级产品包括地面高程和垂直结构参数[23-25]。提取所需地区GEDI光斑内的经纬度、不同百分位树高值(RH)、地面高程、覆盖度、波形参数、信号强度、平均海平面等参数。

1.2.2机载激光雷达数据 机载激光雷达数据(2021年7月)通过规划的航线，使用Riegl·VUX-1无人机激光雷达系统获取平均点密度42 pts·m-2，飞行高度 90 m的数据(飞行区包含设置的28块样地)。获取到机载激光雷达点云数据后，利用点云魔方软件自带的功能处理，获得数字地形模型(Digital Terrain Model，DTM)和数字地表模型(Digital Surface Model，DSM)，通过计算可获得冠层高度模型(Canopy Height Model，CHM)[26-28]。

1.3 研究方法

1.3.1星载激光雷达数据的提取与筛选 根据获取的机载激光雷达数据的经纬度参数，我们在网站进行筛选，并下载所选28块样地星载激光雷达数据[28-29]。提取的参数如表1所示。

表1 提取参数介绍

为了保证获取的星载激光雷达数据更具有，我们除了筛选经纬度，还对其他参数进行了处理。

(1)degrade_flag=0，若该值为“1”，表示指示指向和/定位信息的状态下降，因此删除degrade_flag=1的光斑，保留degrade_flag=0的光斑。

(2)quality_flag=1，质量标志为“1”表示波形满足基于能量、灵敏度、振幅和实时表面跟踪质量的特定标准，可以表示为有效波形。

(3)surface_flag=0，表面标志为“0”表示elev_lowestmode在 DEM 或 MSS 的 300 m范围内，可以表示该参数有效。

(4)sensitivity≥0.9，大于等于0.9的灵敏度阈值表示光斑质量较好，光斑能够穿透毛竹冠层打到地面，因此删除灵敏度小于0.9的光斑。

(5)|elev_lowestmode-digital_elevation_model-mean_sea_surface|<50，如果星载激光雷达数据获取的(地面高程-平均海平面高度)与TanDEM-X的绝对值差值大于50 m，则被认为是受云层影响的错误数据，即删除该光斑。

(6)RH100≤30，通常来说，毛竹不会超过30 m，如果出现了大于30 m的数据则被认为是受云雾影响或非毛竹冠层高度，数据不正确。

(7)landsat_treecover>0，modis_treecover>0，由于毛竹与阔叶树相比冠幅小，激光更易打到地面以造成误差，故删除覆盖度不是大于零的错误值。

1.3.2机载激光雷达数据重采样 机载激光雷达数据较为精准，具有1 m的高分辨率，GEDI光斑直径约为25 m，因此需要对机载激光雷达进行重采样，降低其分辨率至25 m[30-31]。首先，研究使用了ArcMap最邻近法对DTM进行重采样；其次，我们用Matlab软件编写重采样代码，提取CHM不同的百分位树高值(CHM 90～100)；最后，使用ArcMap表面分析工具提取DTM中的坡度、坡向等地形指数数据。

1.3.3星载激光雷达与机载激光雷达数据融合 星载激光雷达与机载激光雷达数据的融合，既能解决机载激光雷达数据获取成本高的问题，又能解决星载激光雷达数据分布不连续的问题。首先运用ArcMap的投影与变换将重采样过的星载激光雷达数据与机载激光雷达数据投影到同一坐标系；其次根据经纬度坐标，将CHM 90～100、坡度和坡向与星载激光雷达数据融合；最后剔除融合数据的重复值与异常值。研究使用70%的大别山毛竹林机载激光雷达数据进行建模，30%的数据进行模型验证。

1.4 模型构建与验证

首先对星载数据RH 90～100与CHM 90～100进行拟合分析，找出拟合程度最优的一种组合，即为机载数据最匹配的百分位树高值，以避免机载激光雷达数据对模型估算产生误差。然后进一步分析星载激光雷达的不同特征参数对毛竹冠层高度精确度的影响，我们共提取5个变量，分别为植被覆盖度、太阳高度角、太阳方位角、坡度与坡向，利用随机森林对5个特征参数进行判别，找出关键影响因子。再利用逐步回归构建的毛竹冠层高度估算模型，根据均方根误差(RMSE)和决定系数(R2)评价建模效果，选出最优模型。最后，利用样地实测的毛竹数据验证模型预估效果。

2 结果与分析

2.1 激光雷达数据融合

我们将网站下载的星载激光雷达数据与无人机获取的机载激光雷达数据在ArcMap中进行融合，计算每个区域的CHM 90～100、高程、坡度和坡向。图2a为研究地冠层高度融合的分布图，由蓝色到黄色依次代表冠层高度由低到高，其中红色的点连成的线是星载激光雷达的激光器发射到地面上的光斑。我们发现，该地毛竹冠层高度在0～25.8 m之间。图2b为星载与机载激光雷达融合的高程分布图，该地地势中间高，四周低，高程差为达到264 m，高程主要在580～640 m之间。图2c为数据融合的坡度分布图，该地中部地势较为平坦，南部、东部及东南部坡度较大。图2d为数据融合的坡向分布图，主要朝向西北方、北方、南方和东南方，其中坡向为西北和东南的占比达到了70%。

图2 研究地多源数据分布图(a)冠层高度分布图 (b)高程分布图 (c)坡度分布图 (d)坡向分布图Fig.2 Distribution maps of multi-source data in the field(a) Canopy height distribution map (b) Elevation distribution map(c) Slope gradient distribution map (d) Slope aspect distribution map

2.2 百分位冠层高度匹配

基于机载激光雷达与星载激光雷达融合的数据，对RH 90～100与CHM 90～100进行相关性分析，结果如表2所示。星载与机载激光雷达冠层高度的数据匹配效果好，R2在0.75以上的达到92.6%。其中，RH 90与CHM 100的R2值最小，为0.63；RH 97与CHM 90的R2值最大，为0.90。对比分析发现，表格数据以RH 97为中心，呈现向两边递减的趋势。包含RH 90、CHM 100的组合匹配程度较低，包含RH 96、RH 97、RH 98的组合匹配程度较高。

2.3 特征参数重要性分析

利用随机森林，分析植被覆盖度、太阳高度角、太阳方位角、坡度与坡向这5个特征变量对毛竹冠层高度精确度的影响。特征变量重要性的大小排序如图3所示，太阳高度角和太阳方位角对星载激光雷达数据精确度影响较大，植被覆盖度以及坡度也是影响星载激光雷达数据精确度的重要原因，坡向的影响相对较小。

2.4 大别山毛竹冠层高度的预估与验证

根据表2，我们发现CHM 90和RH 97拟合程度最高，R2达到了0.90，所以将RH 97以及特征参数分析的植被覆盖度(landsat_treecover)，TAN坡度，坡向，太阳高度角(solar_elevation)，太阳方位角(solar_azimuth)作为自变量，CHM 90作为因变量进行逐步回归分析，建立多个毛竹冠层高度预估模型，算出预估冠层高度与实测数据之间的决定系数(R2)和均方根误差(RMSE)(表3)。结果表明，在加上植被覆盖度(landsat_treecover)这个参数后，R2由0.63提高到0.75，RMSE由2.39 m减少到1.81 m。加上坡度参数后，模型方程的精确度又进一步提高(R2增加了0.08，RMSE减少了0.19 m)。

同时，我们将回归模型预估出的毛竹冠层高度进行验证，如图4所示。与实测数据相比较我们发现，均方根误差(RMSE)达到1.33 m，决定系数(R2)达到0.96。

表2 星载激光雷达数据与机载激光雷达数据匹配结果

表3 毛竹冠层高度拟合方程

图3 特征变量重要性排序Fig.3 Sorting characteristic parameters by their importance

图4 实测数据验证毛竹冠层高度Fig.4 Estimated Moso bamboo canopy heights verified by measured data

3 讨论与结论

激光雷达数据具有较强的观测高度参数提取能力。毛竹是禾本科刚竹属单轴散生型常绿乔木状竹类植物，与阔叶树、针叶树不同，毛竹冠幅小。与冠幅大的阔叶树和针叶树相比，激光雷达发射的激光束更易打到地面[8,13]，如何获取有效数据是研究的主要难点。

基于随机森林算法分析各特征参数对GEDI数据反演冠层高度的影响时，我们发现除与太阳辐射有关的2个参数外，毛竹的覆盖度参数对激光雷达反演毛竹冠层高度的精准度影响非常大[26，32]。影响阔叶树或针叶树等大冠幅树种反演精度的最主要的参数，可能不是覆盖度，反而是坡度等地形因素，这要求分析变量的重要性部分有所不同[33-35]。研究表明在复杂地形条件下需要考虑太阳高度角、太阳方位角、植被覆盖度等因素，以提高毛竹林冠层高度反演模型的估算精度。机载激光雷达空间分辨率为1 m，GEDI星载激光雷达分辨率约为25 m。为了让二者数据进行融合，我们用ArcMap将机载激光雷达数据下采样到25 m[35]。在匹配最优百分位数值时，研究选择了CHM 90和RH 97这一组合，R2和RMSE分别为0.96和1.33 m。Potapov等[36]利用GEDI数据构建全球森林冠层高度模型RMSE为6.64 m；Lang等[37]利用GEDI数据构建全球森林冠层高度模型，该模型产生的全球冠层高度估计值的预期RMSE为2.78 m。因此，利用GEDI数据可以大尺度精准的反演冠层高度。

综上所述，笔者将星载激光雷达数据与机载激光雷达数据进行融合，以毛竹为研究对象对其冠层高度进行反演，在评估特征重要性参数的基础上，通过逐步回归建立冠层高度预估模型。模型中充分考虑太阳方位角、高度角、植被覆盖度、坡度、坡向等因子，与实测数据相比较，模型精确度高(R2=0.96，RMSE=1.33 m)。研究结果为大尺度高精度高分辨率反演毛竹冠层高度提供了技术支撑。