朱琦, 郭华东, 张露,3, 梁栋, 刘栩婷, 万祥星
(1.中国科学院空天信息创新研究院数字地球重点实验室,北京 100094; 2.中国科学院大学,北京 100049;3.三亚中科遥感研究所,三亚 572029; 4.中国林业科学研究院资源信息研究所,北京 100091)
热带森林在生物多样性保护、调节气候变化、防止水土流失、碳汇估算等方面,起到重要的作用[1-2],是自然资源和生态环境的基础。我国热带森林资源丰富,其中海南岛热带天然林是我国连续分布面积最大、物种最丰富的热带林区。2006年,海南岛热带天然林面积为65.90万hm2,主要分布于海南岛的五指山、尖峰岭、霸王岭、吊罗山、黎母山等中部林区[3]。但是由于人类对于热带森林的过渡利用和开发,热带雨林年毁林率逐年增加[4-5]。因此开展热带森林调查,及时准确地获取植被类型覆盖信息,是研究区域内生态环境建设的重点[6]。
遥感对地观测技术的发展,为高效准确地获取实验区域内植被覆盖度提供了有效途径[7-8]。目前,遥感技术可以通过测定森林各项指标特征,协助提取森林覆盖信息、划分森林类型[8]、动态监测森林变化、监测植被生长状态及发育情况等。国内外学者对于如何从遥感影像中自动提取地物信息进行研究[9],也形成了不同的分类体系,取得了丰富的成果。但是目前大多研究都集中在划分森林类型和林种的粗分类,典型热带林种分类精细度依旧较差,无法实现对热带森林资源实时及准确的监测并获取植被覆盖类型信息。
基于多时相遥感数据的影像序列可以反映植被的季相特征以及物候变化规律,有效增加植被的识别精度,目前已经被应用于植被精细提取中[10]。无论从区域尺度的精细分辨到全球尺度的植被覆盖信息提取[11-13],都证明了时序信息的价值。Shimabukuro等[14]利用长达18 a序列的Landsat数据获取了美国犹他州北部因塔特山脉的森林林型分布及变化图,结果较美国森林服务局精度提升60%; Pouliot等[15]提出了一种基于变化检测的时序信息遥感分类方法,获取了加拿大地区2000—2011年土地利用分类图,结果同加拿大国家对于全球土地利用覆盖产品保持一致性和高度准确性。上述研究通过使用多时相遥感数据获取了研究区域内森林林型以及土地覆盖类型的变化。但是目前针对多时相数据的引入相较于单时相数据对于分类精度的具体提升研究较为缺乏,同时对于多时相遥感影像数据分类中的具体策略,如进行影像分类的数量以及时序数据的组合方式等关注较少。
因此本文将多时相信息引入热带森林分类,重点研究如下问题: ①多时相比单时相在热带森林分类中的优势; ②参与分类的多时相数据的数量与分类精度关系; ③参与分类的多时相数据的组合是否影响分类精度。从而提出一种基于时相数据的针对热带天然林的分类方法。由于多时相数据的处理和组合分析对数据数量和处理能力的要求,本文依托Google Earth Engine(GEE)平台,以我国海南岛热带天然林为研究对象,结合长时间序列Landsat8影像,明确遥感影像数量及组合方式对海南岛热带天然林分类精度的影响,提出一种稳定的适合海南岛天然林典型林型的时序分类方法,以期为海南岛热带天然林监测与生态评估等提供有效数据支持。
尖峰岭地区位于海南省西南部乐东黎族自治县和东方市交界(N18°20′~18°57′,E108°41′~109°12′),总面积为640 km2。尖峰岭地区地属低纬度热带岛屿季风气候,常年受西南季风和东南季风的影响,年度雨量分配不均匀,有明显的干湿两季。干季年平均气温最高可达到27~29 ℃,湿季气温较低,仅为16.2~20.6 ℃[16]。该地区自然环境优越,适合植被生长,存在一片我国现有面积较大,保持较完整的热带原始森林[17-18],森林覆盖率达93.18%,类型丰富,包含了山顶苔藓矮林、热带山地常绿阔叶林、热带山地雨林、热带常绿季雨林、热带半落叶季雨林等典型热带植被类型[19-20],是开展热带森林遥感分类研究的理想区域。研究区位置如图1所示。
图1 尖峰岭研究区位置
1.2.1 遥感数据
本文采用GEE平台提供的Landsat8卫星遥感数据开展相关研究。共获取研究区域内2018年度Landsat8 OLI卫星影像29期,行列号分别为124/47和125/47。这些数据是经过辐射校正和几何纠正等处理后的反射率产品[21- 23]。图1中影像是GEE提供的Landsat8影像示意图。
1.2.2 辅助数据
1)海南省土地利用分类图。本文通过中国科学院地理科学研究所与资源研究所网站获取2018年海南省中国土地利用遥感监测数据,空间分辨率为1 km。海南省土地利用类型主要包括耕地、林地、水域及居民地6个一级类型、25个二级类型。将TIFF文件导入GEE中,按编码进行筛选,得到海南省主要林地分布图,后续通过此林地分布图获取研究区域边界,同时结合Google Earth高分辨率影像获取分类样本。
2)海南省数字高程模型数据。为方便后续进行样本选择以及结果验证,本文使用的SRTM(the shuttle Radar topography mission)高程数据,空间分辨率为30 m。观察到尖峰岭林区存在明显的由近海岸地区(平均海拔50 m)至林区腹地海拔逐渐攀升,且东部地区海拔较西部明显偏高的现象。
本文采用的研究技术路线如图2所示,包括数据预处理、森林分类类型确定、特征分析、时相数据分类方法以及验证等步骤。
图2 基于时相变化的海南岛天然林遥感分类技术路线
本文采用的Landsat8数据是经过辐射校正和几何纠正等处理后的反射率产品。此外为了获取1 a内研究区每月一期的遥感数据,对同一月份的多景数据,采用GEE内置的Landsat.simpleComposite()方法进行合成,从而形成该月份的合成数据。合成过程中,该方法使用SimpleLandsatCloudScore()算法标注每个像素受到云的影响程度,然后选择受到云影响最小的像素进行影像合成。经过处理,获得了2018年1—12月每月的合成影像数据。
依据宋永昌[24]在2011年提出的海南岛森林植被分类系统,参考在尖峰岭研究区的实地调查,本文将该地区热带森林分为典型热带雨林、热带季雨林以及常绿阔叶林3种类型,分类级别属于第三级植被型组类型。考虑到研究区域内,常绿落叶林只是零散地分布于常绿阔叶林之中,两者之间并没有明显的区别,所以将常绿阔叶林和常绿落叶阔叶混交林分成一类,统一归类为常绿阔叶林[25]。
基于上述热带森林分类类型,依据考察结果以及高分辨率影像人工判断,选取分类样本。利用这些样本以及获得的光学遥感影像,进一步分析典型热带雨林、热带季雨林以及常绿阔叶林类型的光谱特征、植被指数特征以及时序特征。
2.3.1 光谱特征
光谱特征是最基本的分类依据,基于Landsat8 影像中不同地物类别的典型光谱特征曲线,分析各类别在影像中体现的光谱差异性以及光谱域类别区分性。具体统计结果如图3所示。图3中柱状图代表3类地物在不同波段下的最大反射率与最小反射率的差,折线图表示标准差,两者数值越大代表待分类地物在此波段下区分度较高,有利于后续分类。根据图3可知,3类地物在近红外波段(B5),短波红外1波段(B6)和短波红外2波段(B7)有较大的亮度差及标准差,包含遥感信息丰富。有利于3类地物的区分; 值得注意的是,典型热带雨林在海岸、蓝光、绿光和红光波段都展示出较大差异,分析是由于其季相变化显著,因此波段统计差异较大。
2.3.2 植被指数
植被指数可以有效反映地球表面植被覆盖程度,反映植被的生长状况。其在遥感分类中有助于增强遥感解译以及目标识别能力,目前已经被广泛地应用于土地覆盖类型分类、农作物分析与识别、林地生长与检测。本文中采用的植被指数主要有: 归一化植被指数(normalized difference vegetation index,NDVI)、增强植被指数(enhanced vegetation index,EVI)、归一耕作指数(normalized difference tillage index,NDTI)、陆表水分指数(land surface water index,LSWI)、土壤自适应植被指数(soil-adjusted vegetation index,SAVI)。具体计算公式分别为:
,
(1)
(2)
,
(3)
,
(4)
(5)
式中BLUE,RED,NIR,SWIR1,SWIR2分别代表蓝光、红光、近红外、短波红外1、短波红外2波段影像。图4反映了5类植被指数的统计信息,5种植被指数的标准差、亮度差均明显大于原始光谱信息,其中NDTI与EVI指数包含信息又明显优于其余植被指数,典型热带雨林NDTI指数差异性较大,常绿阔叶林相较于热带季雨林在NDVI指数体现较强差异。
图4 植被指标信息统计
2.3.3 时序特征
由于单时相遥感影像包含信息量有限,无法反映地物类别随时间变化的趋势。同时对于植被分类来说,仅依靠单时相影像无法达到较好区分效果。故将2018年影像按照月份进行分组,研究时序特征对遥感分类精度的影响。
几类典型波段以及植被指数时序统计如图5所示,其中横轴代表月份,纵轴代表该波段反射率。3种地物在光谱特性曲线的变化趋势基本一致,其中常绿阔叶林同其余2类地物区分较为明显,典型热带雨林和热带季雨林在光谱特性中差异并不大。海岸波段、蓝光波段、绿光波段和红光波段在3类地物的变化趋势基本类似,区分度较小。植被指数特征差异较光谱特征明显增大。同时在8—9月,常绿阔叶林与典型热带雨林和热带季雨林在多个波段的反射率有明显变化,此时研究区即将进入旱季,常绿阔叶林开始出现落叶现象,之后进入10—11月,典型热带雨林和热带季雨林反射率均有所下降。基于此,最终选择近红外波段(B5)、2个短波红外波段(B6和B7)以及5类植被指数共计8组变量进行后续分类探究。
(a) 蓝光波段 (b) 短波红外波段
(c) NDVI (d) EVI
图5 典型地物2018年度时相光谱特性曲线
利用前面合成的研究区2018年12期数据,以及热带森林类型的敏感特征集,获取一个年度特征时相信息。考虑到支持向量机(support vector machine,SVM)在处理小样本、非线性、高维数问题的优势,选择SVM作为本研究的分类器。其是由Vapnik于1995年提出,以结构最小原理和VC维理论为基础,通过构造函数子集来达到风险最小化的目标,在大大降低了建模复杂度的同时,提高了分类器学习的能力,具有较高的实用性[26]。
本文设计如下4组实验: ①首先探究单时相影像可以达到的分类精度; ②为了验证多时相数据对于分类精度的影响,探究了不同时相数量下分类精度的变化; ③同时探究了使用筛选后特征与原始特征对于分类精度影响的对比; ④还设计了不同种时相数据的组合方式,如原始单时相分类精度较优或较差的组合,按照季度筛选影像的组合方式等,探究其对最终分类结果的影响。
本研究中分类的精度通过获得混淆矩阵后,计算总体精度(overall accuracy,OA)和Kappa系数进行评估,评估过程中将随机选取70%样本用于训练分类器,剩余30%样本用于验证分类结果。
依据第2.3节波段选择的结果,首先对12期不同时相的Landsat8影像达到的分类精度进行了研究,分类精度如图6所示。
图6 多波段Landsat8影像分类结果
使用多波段进行分类时,在6—8月时期,研究区处于湿季,地区区分度较小,整体分类精度较低(55%); 分类精度最优值在每年的11月份,这与研究区进入干季的时间相同,此次诸如常绿阔叶林和热带季雨林均有落叶现象,而典型热带雨林仍保持常绿状态,同其余2类地物有明显的差异,因此整体分类精度较高。
为了分析多时相数据对分类的影响,将12期不同时相数据依次加入分类器,得到的分类结果精度如图7所示。图中横坐标表示在原有月份基础上,追加当前月份影像数据。
图7 2018年月度时序分类精度变化
根据图7,随着时序影像的加入,OA及Kappa系数均持续提升,加入第二组时序影像时,OA提升最为明显(10%), Kappa系数也有显著的提升; 随着影像的加入,两者增加的趋势逐渐变缓,最终稳定在90%左右。OA最终达到峰值(91.19%),Kappa系数同样到达极值(0.868 5); 最后一组影像加入时,两者均有小幅度的下降。通过对比混淆矩阵,整体而言,仅有小部分的热带季雨林和常绿阔叶林未得到良好区分,出现错分现象,但是时序信息的加入对于遥感影像分类有明显的促进作用。分类精度较单一时相分类精度有明显提升。
为了分析时序影像分类时每次加入影像对于精度的提升,图8表示从第二期影像加入时,每期影像加入时分类精度的提升,随着时序影像序列的加入,OA及Kappa系数的提升逐渐进入平缓期,最终会出现负增长。分析两者变化趋势,在加入第5组数据后,每次加入时序影像对于分类精度的提升都在1~2百分点之内,对分类精度的影响较小。因此当时相数据达到一定数量后,分类精度提升有限,本研究中采用4景不同时相图像可以得到最优的性能效率比。
图8 OA及Kappa系数增加变化
由于参与遥感影像分类的波段间可能会具有较强的相关性,而使用相关性较强的波段进行分类在增加信息冗余的同时也会带来更大的计算量。考虑到本文使用的5种植被指数可能会存在信息重复叠加,因此通过分析5种植被指数的相关性,并对其进行筛选。
为计算5种植被指数的相关程度,需要选取特定月份的数据。本文选取原始单月份分类精度最优的11月影像进行计算,这是考虑到由于11月份影像波段间差异最为显著。得到的5种植被指数相关系数如表1所示。为了获取信息最为丰富的植被指数波段,对每个植被指数同其余植被指数的相关系数进行取绝对值后求和。数值越大表示该指数同其余指数间相关性越强,根据表1可知,NDTI和EVI指数同其余整体相关性最弱,因此将其余3类植被指数去除,仅保留NDTI和EVI,重新进行分类。并与单时相影像分类精度进行对比。对比结果如图9所示。总体而言,同原始分类精度差异不显著,略低于原始结果(平均降低1.5%),但是通过此方法可明显减少冗余波段信息,提高分类效率。
表1 5种植被指数的相关系数
图9 特征波段筛选前后OA对比
为了探究时序影像的加入顺序以及不同组合方式对于分类精度的影响,进一步探究了不同种时间序列组合方式对于分类精度的影响,分组包括: a单时相分类结果从高到低排列; b单时相分类结果从低到高排列; c根据季度不同选取; d选取春季及秋季云层含量由少到多排序。由于前文中得到当时相数据多于4景以后,分类精度提高能力有限,因此所采用的不同排序实验都选择4个时相图像。最终分类精度结果如图10所示,信息统计如表2所示。4组时序分类结果均有较明显提升,a组由于本身分类精度较高,因此提升幅度最小,OA值提升11.78%,Kappa系数增幅25.02%; 其余3种时序组合,提升效果非常明显,OA值提升25%左右,Kappa系数也有50%以上的提升。特别是原始分类精度较差结果以及使用数据更易获取的云层含量较小的分类结果,相较于原始影像,时序结合影像的方式OA值分别提升27.15%和28.83%,达到0.721 5和0.844 6,Kappa系数也分别有66.88%和57.58%的提升。同时,使用季度时序影像可以通过使用包含丰富季相信息的4组影像组合,达到良好的分类效果,OA和Kappa系数的增幅分别为21.77%和47.45%。
(a) 精度较优影像组合 (b) 精度较差影像组合
(c) 季度影像组合 (d) 云量较少影像组合
表2 时序影像分类精度变化
分类结果表明,时序影像通过季度组合、精度较差结果组合以及更易获取的含云量较小影像组合的方式,可以有效提升有限数据源的分类精度,为使用有限数据源进行遥感影像分类提供新的方法。
最终分类结果采用的是12期影像数据达到的最优分类结果,如图11所示。尖峰岭林区总面积640 km2,典型热带雨林、热带季雨林和常绿阔叶林,3类地物占总林区面积比例分别为49.1%,44.0%和6.9%,总体分类结果与预期及先验知识基本吻合,其中,常绿阔叶林主要分布在海拔较低的近海岸地区,随着海拔的攀升,热带季雨林占主导地位,典型热带雨林普遍分布在海拔较高处。这也证明了热带森林空间位置分布同海拔、湿度等自然因素密切相关。
图11 海南岛尖峰岭天然林分类结果
本文主要基于热带森林分类的不精细不准确的现状,结合目前热带地区热带森林类型精确识别与分类的需求,依托GEE遥感大数据处理平台,使用长时间序列遥感影像探究复杂研究地形内,多波段、多时相遥感影像热带森林类型分类方法,同时发展出一套明确的适合海南岛热带天然林分类使用的优化数量及组合方式。本文的主要结论如下:
1)结合Landsat8 OLI光学数据,分析了多波段光谱特征和5类植被指数(NDVI,EVI,NDTI,LSWI,SAVI)等特征对分类精度的影响。结果表明随着具有额外信息的波段数量的增加,分类精度持续提升。结合多波段数据、5类植被指数以及时相信息的引入,均可有效提升分类精度。最终典型热带雨林、热带季雨林、常绿阔叶林3类地物的总体分类精度达到91.19%,总Kappa系数为0.866 2。结果表明尖峰岭地区3类林种空间分布有明显规律,常绿阔叶林主要分布在海拔较低,湿度较低区域; 热带季雨林主要分布在中高海拔地区; 典型热带雨林主要分布在高海拔高湿度地区,且面积最大,是林区内优势林种。
2)重点研究了时序特征对于分类精度的影响,研究结果表明,随着时序信息的不断引入,海南岛天然林的分类精度不断提升; 在不断引入时序信息时,加入第二期时序影像对于分类精度的提升最为明显,可达10%以上; 随着时相数据的不断引入,分类精度逐渐趋于饱和,当加入到第五期时序影像序列后,每次分类精度的提升已经不足1%。
3)对时相信息的组合方式对于分类精度的影响进行探究,分类结果表明,当获取时相数据有限时,通过季度组合以及使用更易获取的含云量较小影像组合的方式,可以有效提升遥感影像天然林分类精度(平均分类精度提升25%左右),尤其是在参与分类的数据单独分类精度较低时,其多时相组合对分类精度的提升更加明显,体现了参与分类数据时相选择的宽泛性。该研究对使用时相遥感数据开展时序遥感影像分类有积极意义。
综上,本文发展的基于GEE多波段多时相的遥感影像热带森林类型分类方法,有效提升了以典型热带雨林、热带季雨林和常绿阔叶林为代表的林型分类精度。该方法可以明确适合海南岛热带天然林分类的优化数量和组合方式,可以有效提升分类精度,改善了传统热带森林类型分类粗糙的现状,对多时相热带森林类型精细分类提供了新的方法和思路。