宋洁，刘学录*

（1.甘肃农业大学资源环境学院，甘肃 兰州730070；2.甘肃农业大学土地利用研究所，甘肃 兰州730070）

森林作为陆地生物圈的主体，因其在全球碳循环中占据重要地位，并能够为人类社会提供如生物多样性保护、气候调节和水资源涵养等诸多生态系统服务而受到广泛的重视［1］。森林植被分布监测，不仅为区域乃至全球的环境变化研究提供了基础信息［2］，也为林业规划和生态保护政策的制定提供了基本的决策依据。在此背景下，探索低成本、高时效、操作方便并具有一定精度保证的森林识别方法对森林的监测与保护具有重要的意义。

近年来，遥感技术已成为大尺度森林识别的首选［3-4］。被动光学遥感数据包括高光谱［5-6］以及多光谱遥感影像［7-8］在目前森林信息获取研究中的应用较为普遍。然而，多光谱和高光谱数据都只能提供物体表面的光谱信息，可能无法分离出光谱相似但结构不同的物体，例如灌木与乔木［9］。激光雷达是一种基于可见光、近红外或短波红外激光束测距的主动遥感技术，它可以通过物体返回脉冲所记录的三维位置提供关于物体的垂直结构信息［10］。也正是基于此特点，激光雷达数据已逐步被单独、或与光谱数据结合使用于土地覆被分类的研究中：Sasaki 等［11］结合机载激光雷达数据与高分辨率影像，采用面向对象的分类方法获得97.5%的土地覆被分类精度以及48.4%的树种分类精度；Zhou 等［12］通过星载激光雷达ICESat/GLAS 数据的波形特征将研究区域分类为树木、建筑物和平坦地面并获得87.2%的总体精度；Liu 等［13］结合星载激光雷达ICESat/GLAS 数据与Landsat 数据，将研究区分类为农田、森林、灌木、水体和不透水表面，总体精度达到91%。在分类过程中，以上研究均发现与单独使用激光雷达数据进行分类相比，将主动式激光雷达数据与被动式光学遥感数据相结合，综合依据地物的垂直结构特征以及光谱特征进行分类，能够对具有相似垂直结构但光谱特征不同的地物进行更详细的区分，从而获得更高的分类精度。

虽然机载激光雷达能够获得更为准确的地物垂直结构信息，但其昂贵的价格以及较大的数据处理量使其较难在大尺度尤其是地形复杂区域应用［14］。与之相比，星载激光雷达如搭载在冰、云和陆地高程卫星（ICESat）上的地学激光测高系统（geoscience laser altimeter system，GLAS）因其全球性的采样策略、较为准确的地理定位以及能够免费获取的数据已成为大尺度森林信息获取的有力工具［14-15］。作为第一颗大脚印激光雷达卫星，ICESat/GLAS 于2003年1月由NASA 发射，到2009年10月终止运行。7年中，GLAS 上的激光器通过40 Hz·s-1的频率发射信号，以直径约为65 m 的椭圆光斑（脚印）照亮地球表面并记录返回的激光能量来获取测高信息，收集森林地区的激光雷达回波数据超过2.5 亿束［16］。作为近年来最常用的激光雷达数据［12］，GLAS 数据被广泛应用于诸如植被高度测算、植被冠层结构研究、地上生物量估计、林木蓄积量估计、建筑物高度测算以及土地覆被分类等诸多研究领域［14］。然而，其在土地覆被分类方向的研究主要集中在地形较为平坦的地区，而针对山地的相关研究较少。这主要是因为GLAS 数据对地面坡度较为敏感，地形起伏会导致其波形失真从而引起对地物垂直结构的误判［17］。另外，在结合光学影像进行分类时，山地地区由于山体掩映造成的地形阴影会导致严重的“同物异谱”和“异物同谱”的现象，也给山地森林范围及类型的识别带来了较大困难［18-19］。而山地森林由于其所受人为干扰较少，往往在区域甚至更大尺度扮演着重要的生态屏障角色，尤其在目前人为土地利用对热带森林的影响日益增大的背景下［20］，对山地森林尤其是北方山地森林范围及类型进行监测显得尤为重要且迫切。

本研究以位处西北干旱地区的祁连山国家公园肃南县段为例，结合GLAS 星载激光雷达数据、Landsat OLI影像、Google Earth 高分辨率影像、DEM 数据以及样地调查数据，综合利用各数据提供的垂直结构特征、光谱特征、季相特征和地形特征，探索基于多源遥感数据的山地森林识别精度提升方法，在将森林与其他土地利用类型进行区分的基础上，进一步识别不同的森林类型（针叶林、阔叶林和针阔混交林），以期为今后山地森林识别提供借鉴。本研究旨在回答以下3 个方面的问题：1）GLAS 数据提供的垂直结构信息能否提升山地森林范围的识别精度；2）不同的地形特征能否提升山地森林类型的识别精度；3）分类时采用不同的波段运算组合其结果是否会有不同。

1 材料与方法

1.1 研究区概况

研究区位于祁连山国家公园肃南县段，位于97°23′-102°21′E，37°37′-39°37′N，总面积约1.26×104km2（图1）。祁连山地处青藏高原东北部，甘肃与青海省交界处，是我国主要山脉之一。该地区分布着大量的森林、灌木、草原和冰川，是整个河西走廊的“水库”，在水土保持和生物多样性保护中发挥着重要作用。

研究区海拔1770～5740 m，坡度0°～58°。属典型的温带大陆性山地气候，昼夜温差大，年平均降水量450～650 mm，年平均蒸发量2134 mm 左右，年平均气温-3～1 ℃，近10年来逐渐升高。植被分布受海拔高度影响从低到高可分为荒漠草原、森林草原、亚高山灌丛草甸、高寒荒漠和冰雪带。现有森林植被主要由青海云杉（Picea crassifolia）、祁连圆柏（Juniperus przewalskii）、山杨（Populus davidiana）和白桦（Betula platyphylla）组成。其他散生树种有油松（Pinus tabuliformis）、蒙古栎（Quercus mongolica）和侧柏（Platycladus orientalis）。青海云杉主要分布在海拔2500～3300 m 的阴坡和半阴坡上，祁连圆柏主要分布在海拔2700～3300 m 的阳坡、半阴坡或半阳坡上，常与高山云杉混交。山杨和白桦是次生林的主要营建树种，山杨主要分布在沟壑、溪流边和海拔2700 m 以下的洪泛区。

1.2 遥感数据

1.2.1 Landsat OLI 数据 搭载在Landsat 8 上的陆地成像仪（operational land imager，OLI），通过一系列如带宽增强、信噪比提高和附加光谱带等技术改进，目前已成为Landsat 系列中最强的卫星［21］，且近年来已被广泛用于区分不同的土地覆被类型甚至不同的森林类型。本研究通过美国地质调查局国家地球资源观测与科学中心网站（http：//GLOVIS.usgs.gov/）免费获取覆盖研究区域的4 景2018年L1T 级Landsat8 OLI 影像，影像基本参数见表1。

表1 Landsat OLI 影像基本参数Table 1 The basic parameters of Landsat OLI images

1.2.2 ICESat/GLAS 数据 ICESat/GLAS 发布了15 个一级和二级数据产品（标签为GLA1 到GLA15），这些产品可在美国国家冰雪数据中心网站（https：//NSIDC.org/data/icesat）免费获取。本研究使用了二级GLA14产品，它提供了有关激光脚印的地理位置、表面高程、波形参数和反射比（如信号的起始和结束位置）等信息。同时，虽然最初设计的GLAS 脚印直径约为65 m，但在实际任务中它们的大小因激光器不同而不同，激光器1、2 和3 的实际脚印直径分别约为110、90 和55 m［10］。因为3 号激光器脚印的椭圆度最小且接近圆形［22］，本研究仅使用来自3 号激光器的数据，以便量化复杂地形对GLAS 波形的影响。最终在研究区获取的GLAS 数据包括L3F（2006年5-6月）、L3G（2006年10-11月）、L3H（2007年3-4月）、L3I（2007年10-11月）、L3J（2008年2-3月）和L3K（2008年10月），共计215 个脚印点（图1）。

图1 研究区位置示意图Fig.1 The location of study area

1.3 外业调查数据

本研究的外业调查于2018年10月进行，主要调查内容包括GLAS 脚印点内植被高度测量和分类精度验证样本实地采样两部分。由于本研究使用的激光雷达GLAS 数据采样时间为2006-2008年，且研究区内起伏的地形会对GLAS 波形产生较大影响从而引起对脚印点内植被高度的误判，因此，首先对部分GLAS 脚印点内实际植被高度进行测量以对GLAS 数据提取的植被高度进行校正和验证。鉴于不同地形坡度对GLAS 波形的影响可能不同［17］，采样前首先将研究区坡度分为5 级：0°～5°、5°～15°、15°～25°、25°～35°、35°～58°，并通过综合考虑GLAS脚印点的可及性等相关因素，分层随机抽样位于植被区域的65 个GLAS 脚印点进行野外测量。其中12 个脚印点位于0°～5°，18 个脚印点位于5°～15°，13 个脚印点位于15°～25°，17 个脚印点位于25°～35°，5 个脚印点位于35°～58°。在GLAS 脚印点使用精度约为3 m 的手持全球定位系统（global positioning system，GPS）接收器进行空间定位，以所选GLAS 足迹中心为圆心建立直径为55 m 的圆形采样区域，采用样方调查法测量脚印点内植被的高度：将每个样地分为4 个圆形子样地，一个位于主样地的中心，另3 个子样地分别位于中心样地360°、240°和120°的方位角，每个子样地的直径为15 m，每个子样地圆心之间的距离为20 m。直接采用手持激光高度计（高度分辨率为0.01 m）测量每个子样地内植被高度，并计算子样地所测植被高度平均值。同时记录脚印点内植被类型、森林类型以及地形特征。脚印点植被平均高度定义为各子样地植被平均高度的平均值。同时在调查过程中拍摄了脚印点周边具有代表性的景观，并记录其坐标以及土地覆被类型、植被类型等信息，以补充土地覆盖分类的精度验证数据。

对用于分类精度验证的实地调查样本，综合研究区各类别土地分布状况、面积比例等因素，采用分层抽样的方法获取典型地类和主要森林类型样本共计72 个。采用GPS 记录样本中心坐标，详细记录其土地覆被类型、植被类型等信息并拍摄照片，同时对森林范围的样地记录其森林类型、树种组成、海拔、坡度、坡向等主要调查因子。确保外业调查所采集样本涵盖研究区域典型地类和所有主要森林类型，并尽量在空间上均匀分布。

1.4 辅助数据

Google Earth 上的Quick Bird 冬季影像被用作分类辅助数据，近年来，Google Earth 高清影像资源一直被用作基于遥感的土地覆盖分类的可靠参考数据源［23］。

本研究采用的DEM 数据是用研究区1∶5 万地形图数据生成，利用ArcGIS 10.2 软件生成研究区的海拔、坡度及坡向分布图，为GLAS 数据的地形校正以及森林类型识别提供辅助信息。此外，研究区2018年度更新的森林资源二类清查数据（包括林分的详细信息以及土地覆被类型信息，以下简称森林资源清查数据）也被用来辅助建立本次分类的训练及验证样本。

1.5 Landsat OLI 数据预处理

对光学遥感影像进行预处理主要是为了保证影像光谱信息真实且空间位置精确［18］，主要包括遥感影像的辐射定标、大气校正、正射校正、地形辐射校正、拼接、配准和裁剪。对4 景Landsat OLI 影像做辐射定标和大气校正，以消除或减少大气分子和气溶胶的散射和吸收对地物反射率的影响，使植被波谱曲线趋于正常。因为L1T文件已经使用地面控制点和数字高程模型进行了几何精校正，所以坐标精度基本能满足中小比例尺的要求。但由于山区地形复杂，为防止影像发生几何畸变，通过DEM 数据对影像进行正射校正，以改正影像倾斜和投影差，消除因山区地形起伏和传感器系统误差等而引起的像点位移。

在复杂的山区地形中，受地形遮蔽的影响，对应遥感图像每一像元的地面所接受的有效光照有很大差别，增加了地物光谱信息提取的难度。因此，为提取地物的真实光谱值，必须先进行地形辐射校正以减弱遥感影像地形效应。本研究选择C 校正模型［24］基于研究区DEM 数据以及影像太阳高度角与太阳方位角信息进行地形辐射校正。随机选择研究区域10%的像元作为样本，用来拟合未校正前地表反射率和局部太阳入射角余弦值之间的线性回归方程，得到斜率和截距，进一步得到经验因子c。采用目视分析［25］和去相关分析［26］评价地形校正的结果，发现此次地形校正在坡度平缓地区校正效果较为明显，校正后阴影区减少，阴坡和阳坡的亮度差异减小。但在高山及丘陵集中区域，在进行地形辐射校正之前，地表反射率与局部太阳入射角余弦值之间的相关系数为0.7149，校正后相关系数为0.5694，说明出现了较为严重的过校正现象［26］。因此就整个研究区域而言，地形校正效果不尽理想，仍需就“同物异谱”以及“同谱异物”的地物光谱值的提取进行进一步的研究。

对经过校正后的4 景Landsat OLI 影像通过直方图匹配算法进行镶嵌，以Google Earth 高清影像作为参考，选取位于道路、河流交叉处以及建筑物点的79 个控制点，采用二次多项式模型进行相对配准，总体均方根误差（total root mean square error）为0.481。最后采用研究区域的矢量范围文件对影像进行裁剪。

1.6 ICESat/GLAS 数据预处理

1.6.1 无效数据过滤 由于GLAS 数据质量易受到云污染、信号饱和以及大气前向散射的影响，因此，依照以下两个标准首先对无效数据进行过滤：1）使用GLA14 产品中记录的检测标志（i_satNdx=0，i_FRir_qaFlag=15）识别无饱和以及无云的GLAS 数据并保留［27］；2）将GLA14 数据中记录的高程（d_elev）与DEM 数据进行比较，由于研究区域内均为低矮建筑物，因此将高差阈值设为大于该地区最大冠层高度（50 m），当GLAS 数据与DEM 数据的高差超过70 m 时，将其视为潜在的低云数据，并将其从分析中删除［28］。

1.6.2 地理定位精度验证 GLAS 数据在研究区的地理定位精度对其能否作为本次研究分类依据至关重要。因此，参照前人研究［29］将GLAS 脚印点的地面高程与DEM 高程值进行比较以评估其定位精度。

首先利用GDAL 开源库，依据GLA14 产品中记录的每个脚印点的地理位置与脚印长短轴及方位角，将GLA14 的HDF5 格式的数据转换为ESRI ShapeFile 矢量格式，并转换为WGS84 坐标系。对每个脚印内的DEM数据所有像素的高程进行加权平均，与GLA14 数据记录的高程值进行比较，根据距脚印中心的距离，使用以下公式对DEM 中每个像素的值进行加权［30］。

式中：W是GLAS 脚印点内每个DEM 像素的权重；a和b分别是GLAS 脚印的半长轴和半短轴，记录在GLA14产品中；（x，y）和（x0，y0）分别是DEM 像素和GLAS 脚印点的中心坐标；x′和y′是DEM 像素沿脚印点长轴和短轴的坐标；α是足迹长轴的方位角，同样记录在GLA14 产品中。在加权平均之前，将脚印点内所有DEM 像素权重之和归一化为1。

共计随机选取30 个GLAS 脚印点对其地理位置精度进行验证，相关系数R2为0.964，均方根误差RMSE（root mean square error）为1.49 m，因此认为本研究区域GLAS 数据的地理定位精度误差较小［31］，能够满足分类的需要。

1.6.3 GLAS 脚印点提取地物高度的校正 在GLA14 产品中，记录了GLAS 原始波形通过高斯滤波器拟合的1 个波形信号开始范围增量（d_SigBegOff）以及最多6 个高斯波峰质心范围增量（d_gpCntRngOff），本研究将波形信号开始范围增量作为地物顶端波形信号，最低两个高斯波峰中较强峰质心范围增量近似表示地表波形信号。对其求差所得的波形幅度（Wext）可以表示从GLAS 中提取的地物高度，这种方法也被广泛应用于从GLAS 数据中提取森林冠层高度［32］。

在平坦地区，Wext可以近似表示实际的地物高度，但在山区，地形起伏会对GLAS 波形形状产生影响，导致从中提取地物高度出现误差。因此，本研究参考前人研究［32-33］，采用一种基于几何模型的地形校正方法，引入平移系数a和b，得到经过地形校正后从GLAS 数据中提取的地物高度为：

式中：d表示GLAS 脚印点直径，本研究中取值为55 m；θ表示脚印点所在的坡度等级的平均坡度值（°）；c表示光速（m·ns-1）；FWHM（full width at half maximum）表示激光脉冲宽度，为6 ns；a，b为平移系数，采用通用全局优化算法（levenberg marquard，LM）通过实地调查GLAS 脚印点的实际测量植被高度拟合得到；HGLAS表示地形校正后从GLAS 数据中提取的地物高度。

由于所使用的GLAS 数据与待分类影像有近10年的时间差距，因此在利用GLAS 数据分类之前首先通过咨询当地林业及土地资源管理部门将10年间新增的植树造林地区、林木砍伐地区、火烧迹地以及道路及房屋建设或拆迁区域的GLAS 脚印点进行剔除，以免造成土地覆被类型的误分。另外，由于GLAS 脚印点提取的地物高度主要用于分类波形结构相似但垂直高度不同的植被（森林、灌木、草地和农田），因此，GLAS 脚印点实地调查主要在植被区域进行。而10年间乔木的生长量会对平移系数a和b的拟合精度产生影响，从而影响从GLAS 数据中提取的除乔木外其他植被高度的精度。通过分析研究区森林资源清查数据可知，研究区乔木林龄范围在9～100年，其中林龄在15年以下的林分仅占总林分的0.29%，且林龄较小的树种多为杨类，区内中龄林及近熟林为区内森林主要龄组，占总林分的97%以上。因此，认为10年间乔木的高度增长量对灌木与乔木的区分影响较小。并在采用GLAS 实地调查数据拟合平移系数之前，根据树龄及树种生长曲线将10年间的树高增加值从实际测量的乔木高度中减去［34］，以减小因时间差距导致的GLAS 数据提取植被垂直高度的误差对分类精度的影响。

1.7 森林范围识别精度比较

为明确GLAS 激光雷达数据提供的地物垂直结构特征对森林范围识别精度的影响，在分类过程中，本研究仅依据光谱特征以及综合依据光谱及垂直结构特征分别建立分类训练样本。根据研究区地表覆盖状况和应用需求，参考《国家森林资源连续清查技术规定》地类划分标准，将项目区土地分类为森林、灌木、草地、农田、裸地、建设用地和水域。

以研究区GLAS 脚印点为分类训练样本综合依据光谱及垂直结构特征识别研究区不同地类。由于本研究使用的GLA14 产品经过统一的滤波和高斯分解拟合处理，其波形形状主要取决于脚印点内物体的垂直结构，因此，首先根据不同波形形状将研究区分为3 类土地覆被类型：平坦表面、建筑物和植被［12］。一般将GLA14 产品仅具有一个峰值（i_numPk=1）的拟合波形定义为平坦表面［28］。将具有两个峰值（i_numPk=2）的GLAS 拟合波形定义为建筑物［12-13］。值得注意的是，有研究证明［12］，两个峰值的GLAS 波形一般为平顶的建筑物，而坡顶的建筑物由于其屋顶表面未在同一水平面会导致GLAS 波形更加复杂。由于研究区内建筑物多为平顶民居或低层的工业建筑，因此本研究将i_numPk=2 的GLAS 数据作为建筑物样本。植被由于其通常有多个目标在同一GLAS 脚印点中或其本身垂直结构（叶片、枝干、冠层等）的复杂性，其波形往往表现出比以上地物更为复杂的特征，因此将具有3 个以上且数量不定的波峰（3≤i_numPk≤6）的拟合波形定义为植被［35］。

根据脚印点光谱特征对以上3 大类土地覆被类型进行进一步细分，其中植被由于其光谱差异较小，但不同植被类型其垂直高度不同，因此依据高度对其进行分类［36］。参考项目区实际情况及前人研究［35］，当HGLAS≥6.0 m，将其分类为森林，当HGLAS为1.5～6.0 m，将其分类为灌木，当HGLAS≤1.5 m，将其归类为草地或农田。

从研究区森林资源清查数据中获取不同类型的土地覆被样本数据，建立仅依据光谱特征的训练样本对研究区影像再次进行分类并将结果与上述分类结果进行对比，以总体精度（计算混淆矩阵）和Kappa 系数作为分类精度评价指标。为避免分类精度是由于增加了分类样本提升的，保持两次分类的训练样本像素数量一致。另外，为验证影像分类时Landsat OLI 影像参与运算的波段对分类结果产生的影响，本研究在两次分类时分别按照OLI 影像分类时识别植被常用的不同波段组合参与运算，并对结果进行比较。

1.8 森林类型识别精度比较

由于不同的森林类型很难依靠GLAS 数据中记录的波形信息对其进行区分［12，36］，而在依据光谱信息区分森林类型时，由于研究区地形复杂，虽然进行了正射校正以及地形辐射校正等预处理，仍难以完全消除地形影响导致的光谱异常、信息畸变和类别混淆。而山地不同坡向往往分布着不同类型的森林，更使得“同物异谱、同谱异物”现象对不同森林类型的识别造成困难。因此，本研究引入季相及地形特征［37-38］对不同类型的森林进行分类，并比较相关特征的加入对分类精度的提升程度。

首先将获取于冬季的Google Earth 高清影像与Landsat OLI 影像进行对比，借助季相特征，根据森林资源清查数据中不同森林类型的面积比例，采集森林类型单一集中连片，面积约近100 m2的林分，随机建立针叶林、阔叶林以及针阔混交林的分类训练样本。然后将通过研究区DEM 数据生成的海拔及坡向分布图与影像叠加并设置一定的透明度，对同一森林类型根据不同海拔及坡向的影像光谱特征进行取样，建立分别位于不同海拔及坡向的针叶林、阔叶林以及针阔混交林的分类训练样本。最后，分别依据随机建立的森林类型训练样本，位于不同海拔的森林类型训练样本、不同坡向的森林类型训练样本以及不同海拔及坡向的森林类型训练样本进行分类以识别不同的森林类型，以总体精度和Kappa 系数作为评价指标对各分类结果的精度进行评价。

2 结果与分析

2.1 GLAS 脚印点提取植被高度校正结果

按照不同的坡度级别，依据调查GLAS 脚印得到的实测植被高度利用LM 通用全局优化算法对方程（5）进行参数拟合，得到各坡度等级的平移系数a和b的最佳估计值分别为：0°～5°，a=1.15，b=2.91；5°～15°，a=0.58，b=1.13；15°～25°，a=0.62，b=0.60；25°～35°，a=0.38，b=0.52；35°～58°，a=0.29，b=0.51。基于Google Earth高分辨率影像随机选取51 个位于森林的GLAS 脚印点（其中6 个脚印点位于0°～5°，14 个脚印点位于5°～15°，11个脚印点位于15°～25°，14 个脚印点位于25°～35°，6 个脚印点位于35°～58°），依据拟合的平移系数计算GLAS 数据提取的森林冠层高度HGLAS。以研究区2018年更新的森林资源清查数据作为独立验证数据，将10年间高度增加值从其记录的森林冠层高度中减去并分别与地形校正前直接从GLAS 数据中提取的波形宽度Wext以及地形校正后GLAS 数据提取的冠层高度HGLAS进行比较，得到各坡度等级地形校正前后的相关系数（R2）和均方根误差（root mean square error，RMSE）（表2）。

表2 不同坡度等级地形校正前后GLAS 提取植被高度的精度比较Table 2 Accuracy comparison of GLAS derived canopy heights before and after topographic correction in different slope gradient classes

地形校正前，随着坡度的增加，GLAS 数据提取树木高度的RMSE 从2.83 m 增加到12.59 m。而地形校正后，不同坡度等级GLAS 数据提取树木高度的RMSE 稳定在1.84～3.91 m，在很大程度上消除了地形坡度对GLAS 波形的影响，从而使得从GLAS 数据中所提取的地物高度可用于下一步的分类研究。

2.2 森林范围识别结果

按照光谱特征对GLAS 脚印点中被分类为平坦地面的道路、裸地和水域进行进一步区分，3 种地物在各波段光谱特征趋势均具有明显差异（图2），因此结合光谱特征对其进行分类，并将道路与建筑物样本合并为建设用地分类依据。

图2 道路、裸地以及水域的平均光谱曲线Fig.2 Average spectral curves of road，bare land and water

可知，OLI 影像上草地和农田在红波段（Band4）和近红外波段（Band5）表现出明显的植被光谱响应“峰谷”特征，且峰值有一定差异（图3），另外，本研究还参照冬季的Google Earth 高清影像，依据项目区耕作特点，对草地和农田的分类进行辅助。

图3 草地和农田的平均光谱曲线Fig.3 Average spectral curves of grassland and farmland

由于GLAS 脚印点内可能存在不止一类土地覆被类型（如混合有建筑物、农田及树木的脚印），因此只有脚印内某地类覆盖其足迹面积50%以上才被视为该地类。按照OLI 影像常用于植被识别的波段组合参与运算，采用目前已被证实分类精度较好的基于支持向量机的监督分类方法分别依据光谱特征以及光谱和垂直结构综合特征对研究区进行分类。基于实地调查数据以及森林资源清查数据，通过分层随机抽样生成在空间上均匀分布的389 个验证样本（共3627 个像元），分别对分类结果进行精度验证（表3）。

表3 不同波段组合及分类方法基于不同特征的分类结果对比Table 3 Comparison of classification results based on different features of different band combination and classification methods

可以看出，不同的光学波段组合运算对分类精度并无影响。同时，相比仅依据光谱特征进行分类，加入垂直结构特征后，分类精度明显提升，总体分类精度提高了10.67%（表3）。为了进一步明确GLAS 垂直结构特征对各地类分类精度的影响，本研究以NIR、Red、Green 波段组合为例比较了基于光谱特征和基于光谱与垂直结构综合特征分类的不同地类的分类精度（表4）。

由表4可知：1）加入垂直结构信息后，森林、灌木及草地的分类精度明显提升，漏分误差（omission error）分别减少了3.38%，24.82%以及31.05%，错分误差（commission error）分别减少了21.47%，14.48%以及5.30%。说明垂直结构信息对光谱特征相似但垂直结构不同的植被的区分具有明显的协助作用。2）建设用地以及水域的制图精度和用户精度均有所提升。这主要是因为研究区建设用地面积通常较小，仅依靠光谱信息分类时，容易与周边地类光谱信息混淆，而GLAS 脚印点垂直结构信息的加入增加了建设用地的位置精度，从而提高了分类精度。垂直结构信息的加入同时消除了由于光线反射或地形阴影造成的水域的错分和漏分现象。3）加入垂直结构信息后，裸地的用户及制图精度均略有降低，这可能是因为项目区的裸地多为砾石地面，同时山地地形的复杂性也增加了GLAS 数据波形的复杂度，因此加入垂直结构信息后反而使分类精度略有降低。

表4 依据不同分类特征的各土地类型分类精度比较Table 4 Comparison of classification accuracy of different land types based on different classification characteristics

为了进一步验证方法的应用精度，在研究区选择面积为3671.14 km2试验区（占总面积的29%），将基于光谱与垂直结构综合特征分类的面积统计结果与研究区2018年度地理国情及影像解译面积统计结果进行对比。

从表5可以看出，基于光谱与垂直结构综合特征分类结果中各地类所占总面积比重与调查及解译结果各地类比重基本一致，7 种土地覆被类型面积平均相对精度为92.09%（表5），说明将Landsat OLI 影像的光谱特征以及经地形校正后的GLAS 数据的垂直结构特征相结合用于山区土地分类能够提高山区森林范围识别的精度。

表5 各地类分类面积统计与对比Table 5 Area statistics and comparative analysis of different classified land types

2.3 森林类型识别结果

依次根据不同森林类型的随机样本，加入海拔信息特征样本、坡向信息特征样本和海拔及坡向信息特征样本分别对所提取的森林按照针叶林、阔叶林及针阔混交林进行森林类型分类。采用基于支持向量机的监督分类方法以NIR、Red、Green 波段组合参与运算，保持每次分类的训练样本数量一致。

根据不同森林类型的面积比例，基于外业调查数据及森林资源清查数据，选取在不同海拔级别（1770～2770 m、2770～3770 m 和3770～4770 m）以及坡向（阴坡、阳坡、半阴坡和半阳坡）上均匀分布的305 个样本点（共2594个像元）对各分类结果进行验证。由于研究区内3770～4770 m 海拔上仅分布少量的针叶林，而4770～5740 m 海拔没有森林分布，因此仅对1770～2770 m 和2770～3770 m 海拔上的森林类型分类精度进行了验证（表6）。

可以看出，地形信息的加入能够显著提升森林类型分类精度，而相比海拔信息，坡向信息的加入对提升分类精度效果更为显著（表6）。这主要是因为地形信息的加入能够尽量全面的考虑到不同地形特征上各森林类型的不同光谱特征，且山体掩映造成的“同物异谱”以及“同谱异物”现象主要由光照区和阴影区造成。另外，虽然森林类型的空间分布格局在区域尺度上主要受到海拔以及纬度的影响，但森林类型分布格局在局部尺度上主要受到诸如微气候、光以及可获取水分等因素的影响，而在受水分条件限制的半干旱和干旱地区森林类型分布受这些因素的影响尤为显著。因此，就本研究区域而言，加入坡向信息比海拔信息更能够提高森林类型的分类精度。

表6 地形信息对森林类型分类精度的影响比较Table 6 Comparative analysis of the effect of terrain information on the accuracy of forest types classification

坡向信息加入后对各坡向森林类型分类精度的提升作用由高到低依次为半阴坡（21.76%）、阴坡（17.78%）、半阳坡（12.67%）以及阳坡（9.52%）。当同时加入坡向和海拔信息时，分类精度比单独使用海拔或坡向信息时均有提升，但并不是两者提升精度的线性加和，这表明海拔和坡向从不同的空间角度上提高了分类精度，同时二者提供的信息可能存在重叠。本研究进一步对不同的森林类型分类精度进行了对比，当依据海拔及坡向特征分类时，针叶林的相对精度最高，主要因为其存在大面积纯林，相对最容易提取，而混交林由于其通常以小片面积分布于山地阴影地区的阔叶林及针叶林过渡地带，所以受地形因素影响较大，相对分类精度最低。将以上两层次的最优精度分类结果合并，得到研究区森林信息及土地覆被（图4）。

图4 研究区森林类型及土地覆被分布Fig.4 Forest types and land cover distribution of study area

3 讨论

本研究综合利用多源遥感数据提供的垂直结构信息、光谱信息、季相和地形信息探索山地森林识别精度提升方法，研究的不确定性主要来自以下几个方面：

由于本研究获取的星载激光雷达GLAS 数据的采样时间与其他数据（样地调查、遥感影像和森林资源清查数据）的采样时间有约10年的差距，因此，在利用样地调查数据拟合地形校正模型的平移系数时，首先将10年间的树高增加值从实测的乔木树高中减去，以避免造成GLAS 数据提取的除乔木外其他植被（灌木、草地、农田）高度的误差。虽然区内中龄林及近熟林占总林分的97%以上，但这一方法仍会对部分幼龄林的识别产生较大影响，给森林范围的识别带来不确定性。另外，在利用GLAS 数据分类之前，为避免10年间新增的植树造林地区、林木砍伐地区、火烧迹地以及道路及房屋建设或拆迁区域的GLAS 脚印点影响分类，对这一部分的脚印点进行剔除，这使得能够用于分类训练样本的GLAS 脚印点减少，也是影响整体分类精度的一个重要因素。

为了便于量化地形坡度对GLAS 波形产生的影响，本研究选择仅来自3 号激光器的GLAS 数据以统一脚印大小。为获取尽可能多的数据以用于分类，导致所获GLAS 数据的季节性差异较大（数据采样时间为2-6月和10-11月）。由于叶片是激光能量作用的主要介质，冬季落叶以及冠层积雪可能引起植被高度估计的误差［39］，从而降低分类精度。另外，为避免云覆盖对土地覆被分类的影响，所下载Landsat OLI 影像的季相差异也会降低不同森林类型的识别精度。

作为第一颗大脚印星载激光雷达数据，GLAS 数据在本研究中的脚印点被认为近似圆形且直径为55 m，其内存在不止一类土地覆被类型的脚印点，往往将其识别为脚印点中的主要土地覆被类型。由此造成的混合像元问题也成为影响分类精度的重要原因之一。

2018年9月，NASA 成功发射了延续早先的ICESat 任务工作的ICESat-2，其上搭载的先进地形激光测高仪系统（advanced topographic laser altimeter system，ATLAS）相比GLAS 具有不同的传感器规格以及不同的采样模式，激光雷达采样强度和空间分辨率均大大提高，对植被结构的监测更加精确［40］。其提供的ATL08 产品，已于2019年春季开始陆续向公众提供2018年10月至今的全球陆地与植被高度数据［40］。虽然本次研究在制定外业调查方案时未能获取到该产品，但在未来的研究中，将利用ATL08 更加密集的采样点与更精确的植被测高数据，与遥感图像光谱信息以及纹理信息相结合，进一步提升山地森林类型乃至树种的识别精度。

4 结论

本研究以祁连山国家公园肃南县段为例，综合利用GLAS 星载激光雷达数据、Landsat OLI 影像、Google Earth 高分辨率影像、DEM 数据以及样地调查数据对山地森林信息进行分层次提取，比较了加入垂直结构特征、季相特征和地形特征后对山地森林范围及类型识别精度的提升作用，结果表明：1）相比仅依据光谱特征进行分类，依据光谱及垂直结构综合特征分类时具有相似光谱特征但不同垂直结构的地物分类精度明显提升，森林范围的识别精度提高。说明GLAS 数据提供的垂直结构信息的加入能够增强不同类别之间的可分性，使土地类型分类及森林范围识别更为有效。2）地形信息的加入能够显著提升森林类型识别精度，而相比海拔信息，坡向信息对提升分类精度效果更为显著。3）多源多时相遥感影像提供的季相特征能够对不同森林类型的识别提供帮助，而不同波段组合虽然对地物增强的效果不同，但其对分类精度几乎没有影响。

另外，虽然GLAS 数据提供的垂直结构信息与光谱信息结合能够提高山区复杂地形条件下中等空间分辨率遥感影像森林范围识别的精度，但由于GLAS 数据对地形较为敏感，在坡度大于10°的情况下，波形失真会导致地物高度的估计精度降低，因此在山区应用GLAS 数据时，必须首先对其进行校正以消除地形对波形产生的影响。同时，由于GLAS 数据能够提供2003-2009年间的地球表面测高数据，因此其对历史影像的分类精度也具有较好的提升作用，能够为分析森林变化情况以及制定自然资源管理政策提供依据。