森林作为陆地生态系统的主体，在减缓气候变化，改善环境及保护生态安全等方面有重要意义。传统的森林类型识别以野外调查为主，成本高、效率低，周期长，时效性和准确难以满足林业发展的需求^[1]。而遥感技术具有成本低、可重复、覆盖范围大等优势，将其应用到森林类型的识别和森林资源信息中具有重要的意义。由于遥感技术可用于地形复杂或难以到达的地区获取森林信息，与传统方法互补，减少实地工作量，所以有必要探索多源遥感数据获取森林类型信息的潜力。

近年来，对于森林类型识别已有大量相关研究，主要包括森林类型分类和树种识别^[2-4]。李梦颖等^[5]利用Landsat-8影像对某天然林区森林类型分类，结果表明结合纹理特征与光谱特征能够提高分类精度。田甜等^[6]从QuickBird数据提取纹理特征、光谱特征、植被指数对福建省将乐林场的优势树种类型提取。Markus等^[7]借助Sentinel-2影像分别用基于对象和基于像元两种方法对德国巴伐利亚州两个林区内7种树种进行分类。Zhu等^[8]使用多时相Landsat TM影像，利用支持向量机方法对阔叶林进行精细分类。Xie等^[9]使用资源3号影像对某林场的4种针叶林和3种落叶阔叶林进行分类，结果表明树冠高度特征可以改善一些树种的分类精度。以上研究大多使用单一数据源的光学影像，而近些年遥感数据源向着全天侯、多分辨率、多传感器等方向快速发展，相关研究认为结合多源遥感数据，充分利用不同数据的波段信息，优势互补，能够更好地识别地物，已有学者将多源遥感数据组合的方法应用到森林资源分类与识别当中^[10-13]。雷光斌等^[14]利用Landsat TM和CCD影像对岷江上游地区森林类型分类，总体精度达到87.05%。刘刚等^[15]结合SPOT5、Landsat TM、CBERS 3种遥感数据对某林场进行森林信息提取，结果表明多源融合数据使各类型分类精度都有所提高。虽然使用多源遥感影像进行森林信息提取的研究较多，但将合成孔径雷达（SAR，Synthetic Aperture Radar）数据引入到森林分类的研究较为缺乏，尤其针对地形复杂、云雾雨雪较多、数据获取困难的山地森林，仅用单一数据源难以满足分类精度需求。SAR卫星不受天气因素的影响，全天候、全天时不间断工作，能够获取云、雨、雾影响下的地表信息，Sentinel-2A影像具有覆盖宽、多分辨率、丰富的光谱特征等优势，本文尝试结合光学影像和SAR数据，探讨多源遥感影像在山地森林信息提取的贡献和作用，促进SAR数据在森林类型分类研究中的应用。

在不同数据源的森林类型识别研究中，众多分类方法的选择也很关键。其中，基于多棵CART决策树构成的随机森林算法在训练样本不足的情况下也能达到较高的精度^[16]，已有部分学者将其应用到森林资源分类研究当中，并取得较好的分类结果^[17-20]。因此本文考虑利用Sentinel-1A、Sentinel-2A，Landsat-8影像及DEM数据，基于随机森林分类算法，对甘肃连城国家级自然保护区内青海云杉和祁连圆柏进行分类试验，并分析不同特征变量的重要性。

甘肃连城国家级自然保护区，处于祁连山山地与黄土高原的重要过渡区域（范围大致覆盖了102°26′~102°55′ E，36°33′~36°48′ N），且位于黄河上游的重要支流——大通河中下游，其地理位置特殊，森林生态系统完整，是开展半干旱地区山地森林研究的典型区域（图1）。研究区内主要地物类型为青海云杉、祁连圆柏、高山灌丛草甸、山地森林草原、裸岩等，特殊的地理位置和地形地貌使保护区境内自然条件复杂多样，为森林资源的多样性提供了多样的生存空间。

研究使用的多光谱数据（表1、表2）是Sentinel-2A和Landsat-8（OLI）影像。Sentinel-2影像具有10 m、20 m、60 m 3种分辨率，包括可见光、近红外、植被红边等共13个波段，研究所用影像从欧空局数据中心（https:∥scihub.Copernicus.eu/dhus/#/home）下载，选择覆盖研究区无云且质量良好（2018年7月17日）的两幅等级为1C的Sentinel-2A影像，该数据是经过正射校正和亚像元级几何精校正的大气表观反射率产品，因此只需进行大气校正。本文使用Sentinel-2A分辨率为10 m和20 m的波段。借助欧洲航空局ESA提供的Sen2cor插件对Sentinel-2A影像进行大气校正，再使用SNAP软件中的最近邻法将分辨率为20 m的6个波段重采样到10 m，并结合4个原始10 m分辨率的波段共得到10个10 m分辨率的波段，转换成ENVI标准格式，在ENVI5.3里进行影像镶嵌和裁剪。

Landsat-8（OLI）影像有9个波段，本实验使用前7个30 m分辨率的波段。从美国地质调查局网站（https://glovis.usgs.gov/）获取一幅研究区范围内无云的Landsat-8影像（2018年5月16日），利用ENVI5.3软件对影像进行辐射定标、大气校正、图像裁剪等处理，将分辨率重采样到10 m。

Sentinel-1卫星是全球环境与安全监测计划（Global Monitoring for Environment and Security）中的地球观测卫星，由A、B两颗极轨卫星组成，载有C波段合成孔径雷达，全天候、全天时不间断运行，无论天气状况如何，都能获取云覆盖下的连续图像。从欧空局数据中心选取干涉测量宽幅模式（IW）的双极化（VV＋VH）的Level-1级GRD格式的Sentinel-1A雷达数据（2018年9月23日）。在SARscape 5.2.1中对获取的研究区Sentinel-1A原始数据进行地理编码、辐射定标，使用Frost滤波器进行滤波以减少斑点噪声，为避免小值影响，将处理结果转换为以dB为单位的后向散射系数图像，如下所示：

其中：N是通过地理编码和辐射定标计算的像素值。利用ENVI5.3对所得的后向散射系数图像进行研究区范围裁剪，并重采样到10 m分辨率。

SRTM DEM(Shuttle Radar Topography Mission DEM)数据是由美国国家航空航天局（NASA，National Aeronautics and Space Administration）和美国国家地理空间情报局（NGIA，National Geospatial-Intelligence Agency）合作采集和发布的数字高程模型，包括SRTM1（30 m）和SRTM3（90 m），该类型数据覆盖范围广、可免费下载，在地学研究领域应用广泛。获取覆盖研究区范围30 m分辨率的SRTM1 DEM，利用ENVI5.3和ArcGIS10.2进行裁剪，提取坡度、坡向，并重采样到10 m分辨率。

使用的样本数据来源于第二次青藏高原综合科学考察研究于2018年9月28日到10月10日野外实地考察，借助无人机与相机拍摄记录青海云杉和祁连圆柏以及其他林地的树木形状特征、分布特征，以此辅助Google Earth高分辨率遥感影像目视解译（图2）。通过以上方式获得保护区内样本474个，用于分类模型建立和精度验证，其中青海云杉80个，祁连圆柏82个，其他林地101个，非林地211个。研究区的林地大部分位于该区中部，而青海云杉和祁连圆柏多分布林区中央、河谷地带，因此如图2所示，青海云杉和祁连圆柏的样本点分布较为集中在整个研究区中间部分。

研究区属于干旱半干旱山地森林区域，其森林结构复杂，青海云杉和祁连圆柏的生境特点明显，两者的空间分布在坡向上有明显差异，青海云杉主要分布在阴坡和半阴坡，而圆柏主要分布在阳坡。图3是青海云杉和祁连圆柏在 Landsat-8和Sentinel-2A影像中的平均光谱特征，可以看出Sentinel-2A中可见光波段光谱差异比Landsat-8明显，在植被红边波段和近红外波段，Sentinel-2A影像对云杉和圆柏的区分能力优于Landsat-8，总体来看，Landsat-8和Sentinel-2A影像中青海云杉与祁连圆柏的光谱差异在近红外波段、短波红外波段及植被红边波段都大于可见光波段。综合考虑云杉和圆柏在空间分布、光谱方面存在的以上差异以及两种树木在树高、冠幅的差别，选取光谱特征、植被指数、地形因子、纹理特征及后向散射系数共计22个特征变量，作为随机森林分类的输入量。

随机森林算法由Breiman于2001年提出，是一种利用多个分类树对数据进行判别与分类的方法，该方法基本原理是利用bootsrap重抽样方法从原始样本选取多个子样本，并逐一进行决策树建模，每棵树单独完成分类后，由投票得出最终输出的分类结果^[14]。本文基于EnMAP2.2.1软件实现随机森林算法的分类，设置提取林地的4种组合方式和优势树种分类的8种组合方式，依次输入包含不同特征变量的组合以及训练样本集，设置决策树数量为100，每个节点的特征数为默认的特征总数的平方根，计算不同组合方式下RFC（Random Forest Classifier）模型，提取林地和优势树种，获取分类精度最高的一组特征组合。

森林优势树种识别是通过分层分类方法获得的，包括第一层的植被提取，第二层的林地提取和第三层的优势树种分类，技术路线如图4。

第一层分类先从预处理的影像提取植被，使用归一化植被指数（NDVI）进行植被提取，再通过比率蓝波段指数去除错误分类为植被的蓝色屋顶。由于研究区河谷地带的居民区存在一些蓝色屋顶的彩钢房，这种材料的屋顶光谱反射特点与植被相似，值与NDVI相近，导致NDVI 提取过程中无法提取纯净的植被，因此需要利用比率蓝波段指数（RBI）去除错误提取为植被的蓝色屋顶。

其中:ρ_nir和ρ_red分别是Sentinel-2A影像的近红外波段和红光波段的反射率。首先使用逐步逼近法进行NDVI的阈值分析，该方法的思想是先在NDVI直方图中确定初始阈值，然后通过目视检查确定植被与影像最佳匹配时的阈值，再利用逐步逼近法对比率蓝波段指数（RBI）进行阈值分析。

其中:ρ_nir和ρ_blue分别是Sentinel-2A影像的近红外波段和蓝光波段的反射率。经多次试验后，设定非植被区移除的NDVI阈值为0.42，蓝色屋顶清除的RBI阈值为2.3，最终提取出纯净的植被层。

第二层林地从植被层提取，通过四次实验评估S2、S1、DEM和L8对林地提取的相对贡献，包括光谱特征、植被指数，纹理特征和地形特征，用这些特征（表4）进行随机森林运算以提取林地，输出林地区域作为第三层优势种分类的掩膜。

第三层分类是区分林地区域的优势种，获取青海云杉和祁连圆柏的空间分布状况。在光谱特征和地形特征基础上添加后向散射特征以改善分类结果，包括归一化植被指数（NDVI）、比值植被指数（RVI）、增强型植被指数（EVI）和差值植被指数（DVI）和后向散射系数等22个特征变量（表5），以林地作为掩膜基于随机森林算法对青海云杉和祁连圆柏分类。

Sentinel-2A数据和地形特征组合提取林地效果最佳，总体精度为98.43%，Kappa系数为0.931 2，提取结果如图5所示。研究区内地物类型多样，植被类型丰富，使用分层提取能提高分类精度。第一层植被提在去除河流、裸岩、建筑用地时也将河谷地区与植被易于混淆的一小部分蓝色厂房屋顶移除。在纯净的植被层中提取第二层林地时去除高山草甸、灌丛、耕地。为了能准确提取山区林地范围，结合地形特征与光谱特征、纹理特征，对训练样本和验证样本比例设定多次实验，获得取了该区域林地范围。林地范围作为优势树种分类的掩膜，其精度影响优势树种的精度，表6为4种组合的林地提取精度，选取精度最高为98.43%的林地提取结果作为下一步对青海云杉和祁连圆柏进行分类的基础，比较多源遥感数据的8种组合的优势种分类精度。

对比分析8种组合的优势树种分类精度（表7）可知，组合1和2都仅使用了光谱特征，但因Sentinel-2A比Landsat 8光谱信息丰富，Sentinel2A具有的3个植被红边波段和两个近红外波段有利于森林类型的识别，且Sentinel-2A比Landsat-8分辨率更高，使得组合2比1获得更高的精度。选取Sentinel-2A数据为基础，逐步添加Landsat-8数据的光谱特征，DEM数据的地形特征，Sentinel-1A数据的后向散射信息进一步计算并分析。组合4结合两种光学影像得到仅使用光谱信息组合中精度最高的分类结果，比使用单一数据源的组合1和2分别提高8.38%和1.55%。与仅使用Sentinel-2A数据相比，当结合DEM数据，总体精度提高2.50％，该结果也证明了青海云杉和祁连圆柏的空间分布受地形因素的影响较大，而事实两者在海拔和坡向方面分布差异明显。组合6和7在包含所有光谱特征和地形特征的基础上分别添加VV和VH后向散射系数，总体精度并没有提高反而稍有降低，且两组的分类精度差异并不明显。当组合8使用光谱、地形、后向散射特征共22个变量时分类精度达到最高，总体精度为92.85%，Kappa系数为0.895 8，比仅使用单一遥感数据源的组合1和2分别提高11.64%和4.81%。

组合8中青海云杉和祁连圆柏的分类精度如表8所示。青海云杉和祁连圆柏是该区域内的优势种，祁连圆柏又是该区域的建群种之一，本文重点关注青海云杉和祁连圆柏的空间分布状况，因此将研究区内油松、落叶松、山杨、白桦等类型归为其他林地，分类结果（图6）可以清晰看出青海云杉和祁连圆柏的整体分布状况和坡向上的分布差异，两者空间分布趋势与祁连山自然保护局所绘制的连城国家级自然保护区植被图中云杉林和圆柏林一致。

随机森林算法在对输入的特征变量进行分类模型建立时还能计算特征重要性，即不同特征变量对最终分类结果的贡献。对分类精度最高的组合8中22个特征变量做归一化特征重要性评分（图7）。

山地森林的森林类型分布与地形因素有很大关系，图7也表明海拔和坡向对保护区优势树种分类的重要性。山区崎岖不平的地表起伏导致影像部分像元的反射率受到影响，阴影问题尤其在春季或冬季的遥感影像中比较严重，选取夏季时间的光学影像，并在分类过程通过加入DEM、坡向等辅助信息，解决阴影对分类结果的影响，避免地形效应^[21-22]。青海云杉和祁连圆柏分类的最高总体精度来自于将VV和VH后向散射系数与Sentinel-2A、Landsat-8和DEM结合，添加VV和VH后向散射系数有助于总体精度小幅度提升，这与Liu等^[23]得到的结果一致。

本文利用Sentinel-2A、Sentinel-1A和Landsat-8共3种多源遥感影像及DEM数据，基于随机森林分类方法，对连城国家级自然保护区内青海云杉和祁连圆柏进行分类，得到以下结论：

（1）与单一影像相比，添加多源遥感影像和DEM数据不论是林地提取还是优势种分类都能明显提高总体精度。Sentinel-2A较Landsat-8分辨率高，波段丰富，在提取林地时添加DEM得到最高精度98.43%，对青海云杉和祁连圆柏分类时，结合地形特征精度有明显提高，基于之上再添加VV和VH后向散射系数精度达到92.85%。

（2）地形特征中对优势种分类贡献较大的特征变量为海拔和坡向，贡献最小的是坡度，这与青海云杉和祁连圆柏的生长特点有关，青海云杉多分布在山谷、阴坡，祁连圆柏基本分布在阳坡。Sentinel-2A的光谱特征S2mean、植被指数DVI、EVI也具有较高的重要性，而由3个植被红边波段计算的植被指数和Landsat-8的植被指数评分接近。虽然结合Sentinel-1A分类总体精度仅有小幅度提高，但VV极化方式的后向散射系数对分类结果仍有较大贡献。

（3）利用多源遥感影像及数字高程数据，提取光谱特征、植被指数、纹理特征和地形特征，设置不同组合方案，对森林类型进行分层提取，将光学影像的光谱信息和SAR卫星数据的微波信息优势互补，获取比单一数据源更好的分类效果。本文方法在森林类型识别方面有显著优势，在环境复杂的山地森林资源调查、优势树种组识别等方面具有较大潜力。

然而本研究在诸多方面还存在不足，在后续研究工作中需完善并进行深刻讨论：

（1）本文仅基于遥感数据多源化，并未考虑多时相，而植被季相节律信息能通过多时相遥感数据获得，因此后续工作加入多时相的Landsat 8或Sentinel-2A影像，有助于分类精度的再度提高。

（2）随机森林方法可估算特征变量重要性，因此能根据特征重要性排序优选特征变量，剔除贡献最低的特征变量后再次分类估算新的特征重要性排序，重复步骤进行多次试验以寻求最高精度的分类结果和特征组合。

（3）本文提出的研究方法在连城保护区取得较好的分类结果，后续研究中可考虑扩展到大区域的祁连山保护区，提取较大范围青海云杉和祁连圆柏的空间分布范围，尝试大面积山地森林类型制图。