作品编号：D959（二等奖）
作品名称：基于LiDAR和资源三号数据的广西高峰林场森林资源三维动态监测
作者单位：福建师范大学地理科学学院
小组成员：王若琦，康晓莹，许媛媛
指导老师：陈耀亮，陆灯盛
中国是一个发展中大国，也是许多国际性公约的签约国，森林面积和蓄积均居世界前列，人工林面积位居世界第一，中国林业的持续发展无疑对世界林业可持续发展具有重要影响。自1999年国家实行天然林保护工程和退耕还林工程后，国有林场的职能也由以采伐为主转变为管护和抚育为主，植被逐年恢复。森林资源动态监测技术有3S技术、抽样技术、数学方法、数据库技术等。遥感技术不受时间、气候的影响，并且能在大范围内进行可重复性监测，已成为森林资源监测的一个重要数据来源。因此，利用多时期的卫星影像对森林资源变化进行分析, 是近几年林业遥感应用的重要领域。机载激光雷达（light detection and ranging，LiDAR）作为一种新型主动遥感技术，既能获取地物的水平结构信息，同时又能获得地物的垂直结构信息，因而可以实现区域地形的高精度模拟。资源三号卫星是集测绘和资源调查于一体的我国第一颗民用高分辨率立体测绘卫星，目前利用资源三号提取数字表面模型（DSM）的精度较高，因此，利用资源三号立体像对数据和机载激光雷达数据对树高的动态变化进行监测，掌握林场的垂直变化规律，较普通遥感数据来说更为有效、精确。同时，资源三号作为高分辨率的多光谱影像，携带的多光谱信息可以用于辅助反映森林植被信息，且由于其能重复对地进行观测，为林地各类资源的动态变化监测提供了更有力的手段。
图1 研究区域
利用遥感技术, 掌握森林资源水平、垂直变化规律, 为森林资源管理提供可作性强的决策依据。本作品所完成的高峰林场各树种的动态变化监测图，可以为林场森林资源管理提供可操作性强的决策依据，易于林场的工作人员掌握林地各个树种的种植规模、面积的变化从而进一步对林场进行管理。关于林地垂直结构上的动态变化信息，本作品中联合ZY-3数据和激光雷达LiDAR数据提取树高变化的方法是优于只单独利用一种数据的：基于遥感航片数据能够精确提取森林立木树高，但目前基于ZY-3航空立体影像，且在无控制点的条件下大范围获取林分平均树高的研究相对少见。这可能主要因为在森林冠层相对较密的林区，被动光学信号较难穿透森林冠层到达地面，因此ZY-3不能获得林下高精度的地形数据。机载LiDAR点云数据可以生成高精度 DEM，弥补ZY-3立体影像在植被覆盖区域提取DEM的不足，且ZY-3数据成本更低，比单独利用LiDAR数据提取树高更加具有实用性。本作品涉及遥感数据处理模块的总体技术流程如下图所示。
图2 作品制作流程图
3.1 数据预处理
辐射定标，将图像的数字量化值(DN)转化为辐射亮度值或者反射率或者表面温度等物理量的处理过程。辐射定标参数一般存放在元数据文件中，ENVI中的通用辐射定标工具(Radiometric Calibration)能自动从元数据文件中读取参数，从而完成辐射定标。大气校正是为了消除大气散射、吸收、反射引起的误差。目前，遥感图像的大气校正方法最常用为FLAAH大气校正。从波谱对比曲线中可以看出，经过FLAASH 校正的影像，基本去除了空气中水汽颗粒等因子的影响，水体波谱曲线也趋于正常。地形校正是指通过各种变换，将所有像元的辐射亮度变换到某一参考平面（通常取水平面），从而消除由于地形起伏而引起的影像辐亮度值的变化，使影像更好地反映地物的光谱特征。本作品利用ENVI地形校正扩展工具—Topographic Correction，从元数据中读取参数，完成地形校正。
图3 数据预处理——ZY-3 MUX 多光谱数据预处理
3.2基于资源三号ZY-3高分辨率多光谱数据的林地机器学习分类及动态变化
使用预处理后的2016、2018两期的ZY-3数据，根据森林二类调查数据和选取ROI作为训练样本，同样使用随机森林分类器对2016、2018年的影像进行计算机分类，进行分类后处理，得出每年的林地分类结果图和各树种的专题图，通过ENVI软件中的Change Detection工具做出2016-2018年的树种动态分布图和地类转换图。
3.3 联合资源三号ZY-3立体像对和LiDAR数据的树高提取
由于在森林冠层相对较密的林区，被动光学信号较难穿透森林冠层到达地面，因此资源三号数据不能获得林下高精度的地形数据，而机载LiDAR点云数据可以生成高精度DEM，且研究地区的DEM在两年内无大改变。本作品利用ZY-3数据的DSM和LiDAR的DEM数据相减，得出2016和2018年的树高分布。
3.4 林地森林资源的三维动态变化分析
垂直结构上，本作品利用ENVI软件分别提取2016年、2018年ZY-3立体像对的正视、后视和正视、前视两种组合的四组DSM数据；分别对比2016年和2018年的正视后视和正视前视两种组合所提取的DSM数据，选取误差较小的组合为之后研究所用，借助LiDAR数据提取的DSM检验精度之后，利用ENVI中的Band Math工具将2018年的DSM和2016年的DSM相减，得到2016-2018年的林场树高变化分布图，之后对树高的变化根据树种进行动态变化分析和成因剖析水平结构上，本作品通过将2016、2018年的林地利用分类图进行加工，利用ENVI中的Change Detection提取地类转化图，从而进一步对土地覆盖类型变化进行动态监测。4.1训练样本优化
训练样本的质量直接决定了分类效果的好坏，选取训练样本时直接在ZY-3影像上进行目视解译存在一定困难，为获得足够数量与质量的训练样本，采用森林二类调查的地类信息，结合软件Google Earth软件，对获取的点状样本进行筛选，在ArcGIS里将xml文件转为shp格式，利用ENVI里的ROI工具将shp文件转为ROI用于分类。4.2.特征变量的提取
本作品从LiDAR生成的DEM影像中选取进行坡度、坡向和高程的计算。将这些因子与遥感变量结合，作为土地覆盖分类的额外变量。再利用ENVI对ZY-3号影像数据的原始4个波段计算从3*3窗口到23*23窗口的8种纹理特征：在基于概率统计的纹理特征提取中选取Mean (平均值)、Variance (方差)指标作为分类特征，在基于二阶概率统计中选取Homogeneity (均匀性)、Contrast (对比度)、Correlation (相关性)、Dissimilarity (异质性)、Second Moment (二阶矩)、Entropy (信息熵)指标作为分类特征。第三类变量是植被指数，植被指数是对地表植被状况的简单、有效和经验的度量，目前已经定义了40多种植被指数:广泛地应用在全球与区域 土地覆盖、植被分类和环境变化，第一性生产力分析，作物和牧草估产、干旱监测等方面。植被指数用作为分类特征，对植被的提取有较好的帮助，对整体监督分类也有较好的帮助。本作品利用ENVI的Band Math提取的植被指数如图：表1 研究所用植被指数
植被指数 | 公式 |
归一化差值植被指数 (NDVI) | (NIR-Red)/(NIR+Red) |
归一化差值水体指数 (NDWI) | (Green-NIR)/(Green+NIR) |
归一化差值红外指数 (NDII) | (NIR-SWIR)/(NIR+SWIR) |
MERIS陆地叶绿素指数(MTCI) | (NIR-VRG1)/(VRG1-Red) |
红边归一化植被指数1(NDRE1) | (VRG2-VRG1)/(VRG2+VRG1) |
红边归一化植被指数2(NDRE2) | (VRG3-VRG1)/(VRG3+VRG1) |
4.3 基于R语言的变量筛选
再利用ArcGIS中的Arc Tools分别将变量值赋予给样本点，后进行对变量的筛选，本作品中使用基于随机森林算法构建的Boruta算法，在R中对变量进行重要性排序，然后通过计算皮尔森（Pearson）相关系数进行相关性分析，最后选择作为分类特征的变量为以下13种：表2 筛选之后的分类特征变量
变量类型 | 具体特征 |
植被指数 | RVI |
光谱变量 | B1-B4 |
地形变量 | 海拔、坡度、坡向 |
纹理变量 | XB4Me, TXB2Hom, TXB3Sec,TXB1Cor, TXB4Cor |
（注释：TX:基于光谱的纹理变量；B1-B4分别是蓝绿红近红外4个波段；Cor:相关性，Ent:信息熵，Hom:同质性，Dis:异质性，Std:标准差，Me:平均值，Sec:二阶矩，Con:对比度；纹理变量后的数字代表窗口大小）
图4 筛选变量Boruta算法
图5 变量重要性排序图
4.3 基于机器学习的随机森林分类算法
最后将这13种变量利用ENVI软件进行Layer Stacking波段合成，作为随机森林分类的输入影像。本研究中通过ENVI软件中的Random forest classification工具实现，随机森林树的数量（Number of Trees）设置为300，特征数量（Number of Features）默认使用“Square Root”方法，其他设置默认，输入影像与训练样本进行分类。由于初始分类时的分类训练对象较为详细，对分类结果按照分类系统进行合并，再利用ENVI做3*3窗口的majority分析，处理后对分类结果进行精度分析。分类结果图如下：
图6 2016年随机森林分类结果图
图7 2016年随机森林分类结果图
4.4基于ZY-3数据DSM提取
利用ENVI软件中的Terrain/DEM Extraction/DEM Extraction Wizard: New工具，分别提取2016年、2018年的正视后视和正视前视两种组合的四组DSM数据；分别对比2016年和2018年的正视后视和正视前视两种组合所提取的DSM数据，选取误差较小的组合为之后研究所用；在提取的DSM和google earth上分别选取50组对应点，利用ARCGIS中的Extract Values to Points工具，在DSM数据上提取50个点的高程值，将google earth上50个点的高程根据其与lidar数据高程值之间的关系进行转换，从而计算其均方根误差。检验所提取的DSM数据的精度，利用和上述类似的方法，比较分析同期ZY-3 DSM与Lidar DSM、不同期ZY-3 DSM。最后做出整体分析。
图8 2016、2018年DSM提取结果
4.5 联合LiDAR和ZY-3立体成像数据的树高提取
由于在森林冠层相对较密的林区，被动光学信号较难穿透森林冠层到达地面，因此资源三号数据不能获得林下高精度的地形数据，而机载LiDAR点云数据可以生成高精度DEM，且研究地区的DEM在两年内无大改变。本作品利用ZY-3数据的DSM和LiDAR的DEM数据相减，得出2016和2018年的树高分布（如图9）。最后在ArcScene中打开裁剪后的ZY-3多光谱影像，并利用得到的树高数据作为Z轴，制作3D树高图，结果如下图所示。
图9 2016年、2018年的树高提取平面图
图10 2016年、2018年的树高提取3D立体图
4.6 森林资源的三维动态监测
4.6.1水平结构的土地覆盖类型变化
打开两个时相的分类结果图。在Toolbox中，打开/Change Detection/Change Detection Statistics，选择前、后时相分类图（Initial State、Final State）。在Define Equivalent Class面板中，如果两个时相的分类图命名规则一致，则会自动将两时相上的类别关联；否则需要在Initial State Class和Final State Class列表中手动选择相对应的类别，点击Ok按钮。
图11 Change detection statistics结果图
4.6.2垂直结构的树高动态变化
利用ENVI中的Band Math工具将2018年DSM和2016年DSM做差，后在ArcGIS中进行制图，结果如图所示：
图12 2016-2018年树高变化图（上为原图，下为密度分割后图）
为进一步对上述结果进行分析，对其进行了快速统计，并对统计结果进行整理，得到如下结果。
图13 2016-2018年DSM差值变化统计图
结合以上图表：2016-2018年的树高变化，主要变化范围是-20 to -10m，发生以上变化的主要原因为——各树种的自然生长、成林的砍伐、以及采伐迹地种植新的树种。-160 to -30 m和30 to 140m 范围的变化值，是由于DSM数据中的错误点，由于所占百分比极低，因此忽略不计。由于在研究区范围内，桉树和八角种植面积比较广，因此，本研究针对桉树和八角进行了重点分析。4.6.3 桉树树高的分析
将小班数据中的桉树种类提取出来，另存为矢量文件，并用该矢量数据对2018DSM-2016DSM的结果进行裁切，根据所得结果对2016-2018年期间的桉树、八角树高变化进行整体分析。
图14 2016-2018年林场桉树树高变化
（1）为进一步分析桉树树高增加和减少的具体情况，在ENVI中利用band math工具，分别将树高增加部分和减少部分提取出来，结果如下图所示。
图15 Band Math提取桉树树高增高、减少部分图
（2）对桉树树高增加和减少部分进行统计分析，得到以下结果。
表3 2016-2018年桉树树高变化像元数统计
表4 2016－2018年桉树树高增加值统计表
统计结果表明，从2016-2018年，桉树砍伐面积占桉树总面积的1/3左右。根据桉树树高变化统计数据，综合考虑其占比极低及实际情况等，将2016-2018年期间树高增加值大于30m的视为存在误差或错误的数据，经过计算，从2016-2018年，桉树树高增加部分树高平均增加值为7.43m。
4.6.4 八角树高分析
（1）将小班数据中的八角提取出来，另存为矢量文件，并用该矢量数据对2018DSM-2016DSM的结果进行裁切，得到2016-2018年期间八角树高变化的数据，结果如下图所示。
图16 2016-2018年林场八角树高变化
（2）为进一步分析八角树高增加和减少的具体情况，在ENVI中利用band math工具，分别将树高增加部分和减少部分提取出来，结果如下图所示。
图17 Band Math提取桉树树高增高、减少部分图
（3）对桉树树高增加和减少部分进行统计分析，得到以下结果。
表5 2016-2018年桉树树高变化像元数统计
从以上图表分析，八角大多未被砍伐，处于正常生长状态，对于树高降低部分，经过与2018年遥感影响对比发现，导致树高降低的主要原因为砍伐。