结合光学与雷达影像多特征的热带典型人工林树种精细识别

黄 翀 张晨晨,2 刘庆生 李 贺 杨晓梅 刘高焕

(1.中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室 北京 100101；2.中国科学院大学 北京 100049)

近一个多世纪以来，随着人类社会经济的快速发展，大规模天然林采伐已成为全球性问题(Paynetal.，2015)；与此同时，出于对木材的大量需求、经济利益驱动以及生态保护需要，人工林大面积种植，现已成为陆地生态系统的重要组成部分。东南亚是全球范围内人工林主要分布和快速扩张地区之一(Paynetal.，2015)，既包括为工业生产提供原材料的速生林种，如桉树(Eucalyptus)，也包括具有较高经济价值的经济作物林种，如橡胶树(Heveabrasiliensis)、油棕(Elaeisguineensis)。据统计，泰国橡胶林和油棕林面积分别从1961年的约40万hm2和950 hm2增至2017年的约300万hm2和76万hm2，增加近800%和80 000%(http:∥www.fao.org)。人工林建设会对原有的地形地貌和生态环境造成一系列消极影响，如土壤肥力下降、水土流失、生物多样性降低等(Locatellietal.，2015)，准确获取人工林类型和空间格局，对森林碳储量估算、生物多样性保护以及可持续森林管理规划等研究至关重要(Georgeetal.，2014)。

遥感在人工林空间分布制图和时间动态监测方面发挥着重要作用，利用机载高光谱数据(Féretetal.，2012)或激光雷达(LiDAR)数据(Shietal.，2018)能够获得高精度的树种分布信息；然而，由于机载高光谱和激光雷达数据获取成本高昂，对于较大空间尺度的人工林识别，卫星遥感数据得到了更广泛应用，许多研究使用光学卫星影像(如MODIS、Landsat、SPOT等)进行人工林监测和制图，特别是针对橡胶林(Lietal.，2012)、油棕林(Broichetal.，2011)、桉树林(Somersetal.，2010)等。考虑到人工林冠层尺寸，较低空间分辨率遥感影像对其识别存在一定困难(Wuetal.，2009)，而中高空间分辨率遥感影像的观测视场与人工林树种或林分冠层尺寸接近，能够有效减少混合像元、提高植被的细节表达(Gessneretal.，2013)，研究表明，基于高空间分辨率遥感影像，结合影像光谱特征、纹理特征及其他辅助特征可有效提高树种分类精度(魏晶昱等，2016；尹凌宇等，2016)。在东南亚热带地区，由于雨季长，持续的云层覆盖严重，很难获得足够的高质量高空间分辨率光学卫星影像。雷达卫星具有较强穿透性，受云雨天气影响很小，可全天时全天候工作，在热带地区土地覆盖制图中得到越来越多重视(Kohetal.，2011)。国内外学者已开始尝试将光学数据(如MODIS、Landsat)与合成孔径雷达(synthetic aperture radar，SAR)数据(如PALSAR)相结合进行人工林制图(Dongetal.，2013；Gutiérrez-Vélezetal.，2013)，根据SAR数据后向散射系数的差异在大类上划分林地和非林地，人工林与其他林地分类更多依赖光谱特征；但由于这2类卫星传感器在空间分辨率和影像获取时相等方面存在较大差异，一定程度上限制了二者协同应用效果。近年来，欧洲新一代哨兵(Sentinel)系列卫星为在高空间分辨率上协同多源数据进行土地覆盖(变化)监测提供了新途径，目前，单独利用Sentinel-1或Sentinel-2数据在土地覆盖(何云等，2019)、树种识别(Immitzeretal.，2016)、作物分类(韩涛等，2018)和扰动监测(Hojas-Gasconetal.，2015)等方面已开展了较为深入的探索，但协同Sentinel-1雷达信息和Sentinel-2光学信息的树种识别研究尚不多见。

鉴于此，本研究以近年来人工林快速扩张的泰国东部地区为例，结合Sentinel-2光谱特征、纹理特征和Sentinel-1雷达后向散射特征，开展光学与雷达影像协同人工林树种分类研究，评价不同数据源和不同特征在树种分类中的作用，获取最优分类策略，以期为热带典型人工林树种精细识别提供新的技术途径。

1 研究区概况与研究方法

1.1 研究区概况

研究区位于泰国东部Chachoengsao(北柳府)、Chon Buri(春武里府)、Chanthaburi(尖竹汶府)和Rayong(罗勇府)交界处(图1)，主要土地覆盖类型包括天然林、人工林、耕地、水体、湿地和不透水面等；属热带季风气候，日照充足，终年炎热，年均气温28 ℃左右，年均降雨量大于1 000 mm。近年来，由于地方经济发展需要，大面积天然林或耕地被人工林取代，以橡胶树、油棕和桉树为代表的人工经济林种植面积不断增加。橡胶树多种植在平原或坡度较缓的山坡上，行宽株密，行距约8 m，株距约3 m；油棕多种植在平原地区，因种植扩张，近年来向高海拔地区的发展趋势明显，采用三角形种植模式，套种间距9 m，以充分利用光照，最大限度提高产量(Basiron，2007)；桉树人工林排列整齐，大小均匀，树冠较小，多种植在较为平坦的地区。

1.2 遥感数据与预处理

Sentinel-2光学数据和Sentinel-1雷达数据从欧空局(https:∥scihub.copernicus.eu/)免费下载。Sentinel-2数据获取时间为2019年1月17日，Sentinel-1数据获取时间为2019年1月12日。利用欧空局提供的SNAP工具对Sentinel-2和Sentinel-1数据进行预处理。

Sentinel-1数据为Level-1级别地距影像(ground range detected，GRD)IW模式双极化数据，空间分辨率5 m×20 m，幅宽250 km，包含VH和VV两种极化模式。对数据进行辐射定标、多视处理、斑点滤波和地形校正，采用Refined Lee散斑滤波器(Zhangetal.，2014)去除斑点噪声的影响，窗口大小设置为7×7。预处理后Sentinel-1影像的像元亮度转换为以dB为单位的后向散射系数：

σ0(dB)=10×lg DN。

(1)

式中：σ0为后向散射系数；DN为像元亮度。

Sentinel-2数据覆盖13个光谱波段，包括10 m空间分辨率的3个可见光波段和1个近红外波段，20 m空间分辨率的3个红边波段、1个近红外波段和2个短波红外波段，以及60 m空间分辨率的海岸/气溶胶、水汽和卷积云波段(本研究未用到)。原始Sentinel-2数据为经辐射校正和几何校正后的L1C标准产品，利用欧空局提供的Sen2Cor插件进行大气校正，获取地表反射率。采用双线性法将20 m空间分辨率的6个波段重采样至10 m空间分辨率，最终得到10个10 m空间分辨率的波段。

1.3 Sentinel-2光谱特征计算

为充分利用Sentinel-2的高光谱分辨率，验证新增红边波段在树种分类中的作用，本研究计算包括15个红边指数在内的33个光谱指数(表1)，加上Sentinel-2的10个原始波段，共计43个光谱指数。

1.4 Sentinel-2纹理特征提取

Haralick(1979)提出的灰度共生矩阵(gray-level co-occurrence matrix，GLCM)是一种公认的比较成熟、有效的纹理特征提取方法，其使用一个空间共生矩阵计算像素值之间的关系，并利用这些值计算矩阵的二阶统计性质。采用灰度共生矩阵提取纹理特征涉及3个重要参数：窗口大小、步长和移动方向。通过多次试验对比分析，本研究设置窗口大小为7×7，步长为1个像元，移动方向取0°、45°、90°和135° 4个方向的平均值，提取Sentinel-2纹理特征。选取4个重要的纹理参数，包括角二阶矩(angular second moment，ASM)、对比度(contrast，CON)、相关性(correlation，COR)和熵(entropy，ENT)参与树种分类。

1.5 Sentinel-1雷达后向散射特征提取

不同极化模式的后向散射系数及相应的差值和比值已被证明对影像分类具有重要作用(Dongetal.，2012)。本研究通过下式计算VV和VH的差值(Diff)和比值(Ratio)，选择VV、VH、Diff和Ratio共4个后向散射特征作为树种分类的Sentinel-1输入变量：

(2)

(3)

1.6 特征组合及分类

采用随机森林(random forest，RF)模型作为树种分类及特征重要性计算和选择的工具。RF是Breiman(2001)提出的一种以决策树为基本分类器的集成学习算法，与其他机器学习算法相比，RF算法能够并行处理高维、海量数据，对训练数据中的噪声或孤立点具有更强的稳定性和鲁棒性(Wurmetal.，2017；Rodriguez-Galianoetal.，2012)，且具有较为突出的特征选择能力，可利用不参与训练的袋外(out-of-bag，OOB)数据估计每个预测变量的重要性，也可利用OOB数据评估OOB score，以确定RF模型的最佳输入特征，减少特征冗余。本研究中决策树数目(ntree)设为800，每个节点的特征数(mtry)由下式计算：

(4)

式中：p为用于分类的特征数。

由于Sentinel-2包含丰富的原始光谱波段及其衍生出的众多光谱指数特征，全部特征参与分类可能造成信息冗余，导致分类精度降低、分类速度下降，因此，本研究在RF模型中采用平均不纯度减少算法(Breiman，2001)对光谱波段和光谱指数进行重要性评估，根据OOB score确定最优光谱特征组合。

为评估Sentinel-2光谱特征、纹理特征和Sentinel-1雷达后向散射特征对人工林树种识别的能力和贡献，本研究考虑4种特征组合(表2)作为RF分类器的输入特征，比较不同类型遥感数据及其特征在分类精度提升上的表现。

表2 分类特征组合Tab.2 Feature combinations for classification

1.7 精度评价

2018年9月，对研究区土地覆盖类型进行实地调查，获得251个实地样本点。通过Google Earth的高分辨率卫星影像(http:∥earth.google.com/)和Global Croplands的实地图像(https:∥croplands.org)对样本点进行补充，尽可能使样本点均匀分布在研究区内。最终共选取1 397个验证样本，其中包括203个天然林样本、207个橡胶林样本、209个油棕林样本、131个桉树林样本、251个耕地样本、143个水体样本、50个湿地样本和203个建筑用地样本。采用混淆矩阵对4种特征组合的分类精度进行评价，评价指标包括生产者精度(producer’s accuracy，PA)、用户精度(user’s accuracy，UA)、总体精度(overall accuracy，OA)、Kappa系数以及PA和UA的调和平均值(F1)。F1的取值范围为[0，1]，值越大，分类效果越好，反之则越差。F1计算公式(Baumannetal.，2012)如下：

(5)

2 结果与分析

2.1 不同光谱特征在人工林树种分类中的重要性

采用平均不纯度减少算法计算得到Sentinel-2影像15个红边指数、18个光谱指数和10个原始波段的特征重要性如图2所示。由绿波段B3和近红外波段B8A计算得到的NDWI2的重要性最高，为5.02%；其次为绿波段B3和短波红外波段B11，重要性分别为4.45%和4.29%；红边波段B5和红边指数NDVIre1、NDre2和NDVIre2的重要性较高，分别为2.89%、3.97%、3.67%和2.68%，说明红边波段对树种分类具有较高价值；而原始波段B8和B6的重要性最低，不足1%。

图2 Sentinel-2光谱特征重要性排序Fig.2 Feature importance of spectral features for Sentinel-2 imagery

43个光谱特征组合模型的OOB score分析如图3所示。特征数从1增至17，OOB score逐渐增加，特征数为17时达最高值0.947 2，此后有轻微降低趋势，因此本研究选择重要性排名前17的光谱特征参与树种分类。统计重要性排名前17的光谱特征中包含Sentinel-2各波段的特征数，如图4所示。在前17个光谱特征中，包含近红外波段B8A的特征最多(8个)，其次为包含可见光波段B3(7个)、B4(6个)的特征。此外，包含红边波段B5的特征也较多，进一步说明新增红边波段在树种分类中具有较高价值。

图3 Sentinel-2不同光谱特征组合模型的OOB scoreFig.3 OOB score of different feature combination models for Sentinel-2 imagery

图4 重要特征波段统计Fig.4 Statistics of selection for each Sentinal-2 band

2.2 不同人工林树种在纹理特征和后向散射特征上的可分离性

分别提取天然林、桉树林、油棕林和橡胶林在纹理特征和后向散射特征上的特征值(图5)，4种林地均具有一定可分离性。从Sentinel-2纹理特征(图5a)看，天然林在ASM特征上明显小于人工林，在CON和ENT特征上明显大于人工林，可分离性较好；在COR特征上与人工林区分度较小。对于3种人工林，油棕林在CON和ENT特征上值最大，在ASM特征上值最小，COR特征值介于桉树林和橡胶林之间；桉树林在ASM、COR特征上值最大，在ENT特征上值最小；橡胶林在CON和COR特征上值最小，ASM和ENT特征值介于桉树林和油棕林之间。这说明桉树林、油棕林和橡胶林在4种纹理特征上具有一定差异性。从Sentinel-1雷达后向散射特征(图5b)看，4种林地均具有一定可分离性，在Ratio特征上的区分度大于Diff特征。总的来说，Sentinel-2纹理特征和Sentinel-1雷达后向散射特征可作为树种分类的有效特征。

图5 天然林和不同人工林在纹理特征(a)和后向散射特征(b)上的特征值(均值±标准差)Fig.5 Statistics of textural(a)and backscattering features(b)of natural forest and different plantations(mean±SD)图5a的主坐标轴表示ASM、CON和COR特征值，次坐标轴表示ENT特征值；图5b的主坐标轴表示VV、VH和Diff特征值，次坐标轴表示Ratio特征值。The primary axis of Fig.5a represents the values of ASM,CON and COR features,and the secondary axis represents the values of ENT feature;the primary axis of Fig.5b represents the values of VV,VH and Diff features,and the secondary axis represents the values of Ratio feature.

2.3 不同特征组合树种分类结果比较

4种特征组合的分类结果如图6所示。仅利用Sentinel-2光谱特征(S)分类(图6a)，“椒盐现象”非常明显，加入纹理特征(S+T)、后向散射特征(S+SAR)以及结合纹理特征和后向散射特征(S+T+SAR)(图6b-d)后，各地类分类图斑的破碎度均有所降低，“椒盐现象”得到改善。3种人工林的分类图斑较为规整，其中，橡胶林分布最连续，面积最大，主要分布在平原或坡度较缓的山坡上；油棕林多与橡胶林镶嵌分布，地块面积较大；桉树林地块面积最小，主要分布在较为平坦的地区，多与耕地交错分布。

图6 不同特征组合的分类结果Fig.6 Classification results of different feature combinations

不同特征组合分类结果的混淆矩阵如图7所示，精度验证结果如图8所示。对比4种特征组合的分类精度，仅利用光谱特征(S)分类，3种人工林的分类精度均较低，不同人工林之间存在不同程度混淆，且易与天然林和耕地混淆；橡胶林和油棕林的生产者精度不足0.70，桉树林的生产者精度仅0.53，F1仅0.61，总体分类精度为0.75，Kappa系数为0.71。加入纹理特征组合(S+T)分类，天然林与3种人工林之间混淆减少，人工林误分为耕地比例减小，桉树林、油棕林和橡胶林的生产者精度分别提高0.16、0.15和0.08；3种人工林的整体分类精度均有所提升，桉树林、油棕林和橡胶林的F1分别提高0.16、0.11和0.04，总体分类精度提高至0.80，Kappa系数提高至0.77。加入后向散射特征组合(S+SAR)分类，桉树林、油棕林和橡胶林相比仅利用光谱特征(S)的分类精度同样有所提升，F1分别提高0.18、0.15和0.08，总体分类精度提高至0.84，Kappa系数提高至0.81。结合Sentinel-2光谱特征、纹理特征和Sentinel-1雷达后向散射特征(S+T+SAR)分类，总体分类精度和Kappa系数达到最高值，分别为0.85和0.83；其他地类误分为人工林的数量最少，3种人工林均达到最高用户精度；桉树林、油棕林和橡胶林的F1均大于0.80，相比仅利用光谱特征(S)分类分别提高0.19、0.15和0.10。

图7 不同特征组合分类结果的混淆矩阵Fig.7 Confusion matrix of classification results of different feature combinations

图8 不同特征组合分类精度对比Fig.8 Accuracies of classification results of different feature combinations

综上可知，4种特征组合的分类结果中，S+T+SAR的总体分类精度和Kappa系数最高，其次为S+SAR和S+T，S最低。桉树林、油棕林、橡胶林、耕地、水体和建筑用地在S+T+SAR中均取得最好分类结果，天然林和湿地在S+SAR中分类精度最高。

3 讨论

本研究基于Sentinel-2光学影像丰富的光谱特征、纹理特征并结合Sentinel-1雷达后向散射特征，开展热带典型人工林树种精细识别，结果发现，Sentinel-2的短波红外波段B11在树种分类中重要性较高，与Immitzer等(2016)和Sothe等(2017)基于Sentinel-2A进行树种分类时得到的结论一致。在随机森林算法选择的17个光谱特征中，Sentinel-2的红边波段和红边指数所占比例较大，许多学者(Immitzeretal.，2016；Sotheetal.，2017；Mngadietal.，2019)同样强调了红边波段对林地分类的重要性。这是由于植被在红边波段敏感度高(Peerbhayetal.，2014)，尤其是红边波段与叶绿素含量、生物量和冠层结构等反映植被叶片性质的信息密切相关(Dubeetal.，2014；Ramoeloetal.，2015)，因此对不同树种的区分具有重要意义。此外，Sentinel-2光学影像的红波段和近红外波段在树种分类中同样表现出较高重要性，这可能是因为红波段叶绿素吸收能力强，近红外波段吸水能力强，从而提高了不同林地之间及与周围地物的光谱可分性。Nooni等(2014)基于Landsat-7提取油棕林时也发现在红波段和近红外波段油棕的可分离性最高。

仅利用Sentinel-2光谱特征分类，橡胶林、油棕林和桉树林3种人工林的识别精度均较低。由于光谱特征相似，人工林与天然林、耕地混淆较多，不同人工林之间也存在不同程度混淆。Dian等(2015)研究表明，结合光谱特性和纹理特征可有效解决光谱混淆造成的分类误差，提高树种分类精度。本研究在光谱特征基础上加入纹理特征(S+T)后，3种人工林的分类精度均有不同程度提高。不同人工林具有不同的株行距和不同形状、大小的树冠，在遥感影像上会形成特有的纹理，而天然林在树冠大小、树高和密度等方面具有更大变异性，纹理特征加入可有效减少“椒盐现象”，提高斑块的完整性和分类精度。如仅利用光谱特征分类，缓坡上的橡胶林因海拔、坡向、坡度等因素影响呈现不同的光谱特征，易被错分为耕地和桉树林；加入纹理特征后，可有效降低山体效应带来的光谱差异，减少橡胶林与桉树林间的误分，橡胶林斑块也更为完整。

加入后向散射特征(S+SAR)后，3种人工林的分类精度同样有所提升，说明Sentinel-1雷达后向散射特征可有效提高人工林树种间的区分度。Torbick等(2016)在提取缅甸和加里曼丹的橡胶林和油棕林时发现，Senitnel-1雷达后向散射特征对橡胶林和油棕林提取具有重要作用。雷达数据基于地物后向散射特征获得不同于光学遥感的影像，且雷达信号具有一定穿透力，可获取植被表面信息和地表粗糙度等额外信息，尤其是合成孔径雷达(SAR)对橡胶林、油棕林和桉树林在生物量、密度和垂直分层等森林结构信息上的差异敏感性较高，使得SAR数据在树种识别中具有独特优势。此外，不同树种冠层含水量也可通过影响介电常数影响雷达信号(Oonetal.，2018)，能够提高不同人工林树种间的区分度。因此，SAR数据可作为光学影像的有益补充，对光谱特征接近的人工林和作物进行区分。研究区耕地多种植水稻(Oryzasativa)，淹水期长，SAR后向散射特征对土壤湿度和淹没程度十分敏感(Chatziantoniouetal.，2017)，可以弥补光学遥感区分人工林和水稻的不足，提高人工林分类精度。

近年来研究发现，集成多种观测类型和数据模式的多传感器与数据融合技术有助于提高森林制图的细节和精度(Hydeetal.，2006；Descleeetal.，2013)。本研究中，仅利用Sentinel-2光谱特征分类精度较低，加入光学纹理特征和雷达后向散射特征(S+T+SAR)后，3种人工林的分类精度均达到最高，这说明单一传感器或特征进行人工林树种分类仍有很大局限，尤其在生态系统复杂、植被覆盖度高的东南亚地区；而多数据源、多特征结合可充分利用地物在不同传感器上感知的信息，增加不同人工林之间以及人工林与其他地物之间的区分度。本研究结合Sentinel-2光学数据和Sentinel-1雷达数据的人工林树种分类精度明显优于以往仅利用Sentinel-2数据的树种分类精度(Immitzeretal.，2016)，可见多传感器、多特征结合应用的重要性。但需要说明的是，本研究Sentinel-1数据只提取了后向散射特征，实际上由于雷达对植被冠层的穿透能力和对土壤含水量的敏感性，高空间分辨率Sentinel-1数据同样含有丰富的纹理信息。此外，Sentinel-1、Sentinel-2数据均具有较高的重访周期，如何结合其多样的纹理特征和时序特征开展精细化人工林制图值得进一步研究。

4 结论

本研究以近年来人工林快速扩张的泰国东部地区为例，结合Sentinel-2光谱特征、纹理特征和Sentinel-1雷达后向散射特征，采用随机森林算法对Sentinel-2不同光谱特征在树种分类中的重要性进行评估，并分别针对光谱特征、纹理特征和后向散射特征的不同组合进行3种人工林树种识别研究，最终实现该地区3种人工林树种的精细提取，得到如下结论：

1)在Sentinel-2光学影像的43个光谱特征中，选择重要性排名前17的特征参与分类时，OOB score达最高值0.947 2；Sentinel-2的蓝波段、红边和近红外波段及其相应的植被指数在树种识别中重要性较高。

2)桉树林、油棕林、橡胶林和天然林在Sentinel-2影像的ASM、CON、COR和ENT纹理特征和Sentinel-1影像的VV、VH、Diff和Ratio后向散射特征上均具有一定差异性，可作为树种分类的有效特征。

3)仅利用Sentinel-2光谱特征对不同人工林树种的区分度有限，桉树林、油棕林和橡胶林的F1分别为0.61、0.74和0.70。结合光谱特征、纹理特征和后向散射特征，3种人工林的分类精度达到最高，F1均大于0.80，同时其他地物的分类精度也有较大提高，总体分类精度和Kappa系数均达到最高值，分别为0.85和0.83。