基于Sentinel-2和Landsat8 OLI数据融合的土地利用分类研究

赵亚杰,王立辉,孔祥兵,阴海明,池 泓,黄进良

(1.水利部黄土高原水土流失过程与控制重点实验室,河南 郑州 450003；2.中国科学院测量与地球物理研究所环境与灾害监测评估湖北省重点实验室，湖北 武汉 430077；3.中国科学院大学，北京 100049;4.黄河水利科学研究院，河南 郑州 450003)

土地利用变化是全球变化的主要驱动力之一，是全球变化与可持续发展研究中的热点[1]，通过土地利用分类，既可以了解各种土地覆被类型的基本属性，进行农作物种植面积提取、产量预估[2]，还可以认识土地利用的区域结构和分布特点，为进一步研究地域差异奠定基础[3].遥感是土地利用与土地覆被变化研究技术体系中的主要组成部分，目前，具有多光谱与多时相特征的遥感技术已成为获取土地覆被动态监测数据和土地覆被分类的最有效工具[4].常见的遥感影像有空间分辨率较低的MODIS数据，也有中空间分辨率的Landsat系列，以及高空间分辨率的高分一号、QuickBird等[5].其中，Landsat8 OLI影像获取的多光谱影像包含15 m分辨率的全色波段和30 m分辨率的8波段影像，Sentinel-2影像则包括10 m空间分辨率的4波段影像、20 m分辨率的6波段影像和60 m分辨率的3波段影像.不同的影像在土地利用研究中已被广泛应用，李静等[6]利用Landsat影像对大棚菜地信息进行提取并取得较好的地物识别效果，李晓慧等[7]利用Landsat8 OLI时序遥感数据结合地物细微光谱特征对农作物进行分类，并获得较高的总体分类精度，张卫春等[8]基于Sentinel-2影像和地形数据在低山丘陵区采用随机森林方法获得的土地利用分类结果优于传统分类方式所获得结果，为地形复杂区域的土地利用分类提供了参考，杨斌等[9]同样使用Sentinel-2A卫星数据提取了四川黑水河流域的干旱河谷，结果表明Sentinel-2A数据丰富的光谱特征能够较好地满足国土资源普查和生态环境变迁分析.但对于Sentinel-2和Landsat数据融合及其在土地利用中的应用等方面，国内的研究很少，国外学者如Gemine et al[10]也只是针对多种影像融合模型进行了质量评价.

鉴于此，本研究以山东省平邑县为研究区，采用Sentinel-2和Landsat8 OLI遥感数据作为研究区数据源，探讨在基于亮度的平滑滤波(smoothing filter-based intensity modulation,SFIM)、高通滤波(high-pass filtering,HPF)、小波变换(a trous wavelet transform,ATWT) 等3种不同融合模型下，Sentinel-2的10 m高分辨率波段与Landsat8 OLI 30 m波段数据的融合效果，并以最佳的融合波段与融合算法，获取整个研究区的融合影像，采用随机森林机器学习算法分别对Sentinel-2原始影像、融合获得影像和Landsat8 OLI影像进行分类，从分类精度上进一步探讨该融合方法在土地利用分类研究中的适用性.

1 研究区概况及数据源

1.1 研究区概况

图1 研究区2017年Landsat8真彩色合成影像Fig.1 Landsat 8 true color composite diagram of the study area in 2017

研究区平邑县位于山东省中南部、临沂市西陲(117°25′—117°56′E,35°07′—35°43′N)，东西最大横距47.35 km，南北最大纵距66.75 km，总面积约1 824.97 km2.研究区属季风区域大陆性气候，具有冬季寒冷、夏季炎热、光照充足、无霜期长的特点，同时具有明显的山区特征，地势南北高，中间低，略向东南倾斜，中部谷地、陆地土质肥沃，田地分布较多，境内土地覆被类型多样.研究区影像如图1所示.

1.2 数据来源及预处理

遥感数据包括Sentinel-2 L1C级别数据的蓝、绿、红、近红外4个10 m空间分辨率波段，Landsat8 OLI数据的对应4个波段，成像时间均为2017年9月14日.针对Landsat8数据，利用ENVI进行辐射定标和大气校正等预处理.Sentinel-2数据则利用哨兵数据应用平台(sentinel application platform,SNAP)中的sen2cor模块进行大气校正、产品格式转换等预处理.对两种预处理后的影像，再利用ENVI软件裁剪出图1红色区域的子研究区，以便后续研究的应用.

2 研究方法

首先在研究区内随机选取面积约为5.76 km2的子研究区，以SFIM、HPF、ATWT 3种不同融合模型为基础方法，利用Sentinel-2的4个10 m高分辨率波段分别与Landsat8 OLI 30 m波段数据进行融合，从而得到最佳的融合波段与融合算法，在此基础上获取整个平邑县研究区的融合影像，根据2017年8月份平邑县土地调查数据和谷歌地球高分辨率影像获取的样本，参考基于遥感技术的全国生态系统分类体系[11]，并采用随机森林机器学习算法分别对Sentinel-2原始影像、融合获得影像和Landsat8 OLI影像进行分类，最后进行分类精度评定.

2.1 数据融合算法

目前，已经有多种图像融合算法在现实中广泛应用，本研究选取SFIM，HPF和ATWT 3种算法作为融合基础.SFIM算法是一种基于空间域滤波变换的融合算法[12]，该算法简单易用,融合影像的光谱几乎不失真,具有较高的光谱保真度,但在高频空间信息融合上略逊色于Brovey变换方法[13].SFIM变换的融合算法公式：

式中,BS为用该算法生成的融合影像；i为波段值，j和k分别为影像的行、列值；Bl为低分辨率影像，在此代表Landsat的4个波段(30 m)；Bh为高分辨率影像，即为Sentinel-2波段(10 m)；Bm为模拟低分辨率影像，它可以通过对10 m波段进行低通均值滤波处理来获得.

HPF是一种基于高通滤波的影像融合算法[14]，HPF融合算法能够较好地提取全色影像细节信息，并把全色影像的细节信息直接叠加到多光谱影像上[15].融合后影像的边界清晰度增强，可更好地保留光谱特征[16].HPF算法的定义公式：

HPK=Wa×PANKH+Wb×MULKL

式中，Wa、Wb为权，且满足Wa+Wb=1，PANKH为Sentinel-2高空间分辨率影像波段进行高通滤波的结果，MULKL为Landsat8影像进行低通滤波的结果，HPK为锐化后的输出影像.

ATWT是基于小波变化的图像融合算法[17]，该算法利用小波变换将源图像分解到不同频段的特征域中，然后在不同的特征域中采用不同的融合规则和融合算子进行图像融合处理，从而形成新的小波金字塔结构，最后利用小波反变换技术获得融合图像，所以小波变化算法能够对图像细节及边缘信息进行较好的体现[18].

2.2 随机森林分类

随机森林是属于集成学习的一种组合分类算法，集成学习的核心思想就是将若干个弱(基)分类器组合起来，得到一个分类性能显著优越的强分类器.如果各弱分类器之前没有强依赖关系、可并行生成，就可以使用随机森林算法.随机森林利用自主抽样法从原数据集中有放回地抽取多个样本，对抽取的样本先用弱分类器—决策树进行训练，然后把这些决策树组合在一起，最终的分类或预测结果由多颗树分类器通过投票得出[19].随机森林算法较强的抗干扰能力，表现性能的优异，相较于人工神经网络、决策树、支持向量机等机器学习方法，在分类精度和分类效率方面更具优势[20],因此，越来越多的学者在进行遥感影像分类研究时采用随机森林机器学习算法[21-22].

3 结果分析

3.1 融合结果分析

首先在子研究区中，利用SFIM、HPF和ATWT 3种算法，分别以Sentinel-2影像的蓝、绿、红、近红外波段作为高分辨率影像参与到融合中，采用不同算法时各波段的融合效果如图2、图3、图4所示.从图中的视觉效果来看，融合影像能够将Sentinel-2影像波段的空间信息和Landsat8多光谱影像的色彩信息集成在一起，相比于原始的多光谱影像，融合影像在保留其波段信息前提下拥有了Sentinel-2影像的细节信息，同时在色彩方面，融合结果与Landsat8 OLI影像的色彩差距并不明显.但在空间细节方面，利用SFIM和HPF算法获得的融合影像中，道路及人工表面等地类的纹理模糊感较强，边界信息模糊，而图4的融合影像表明ATWT算法对细节的处理明显强于前两种算法，融合的结果也更加精确.以Sentinel-2不同波段作为高分辨率融合影像时，融合效果也不同，可以明显注意到，不论采用哪种算法，近红外波段的融合影像都丢失了很多细节信息，目视效果较差.蓝、绿、红、近红外4个波段的融合效果采用相关系数(correlation coefficient,CC)、光谱角制图(spectral angle mapper,SAM)和相对整体维数综合误差(relative global-dimensional synthesis error,ERGAS)3个指标作定量比较分析.相关系数能够反映两幅影像间灰度信息的相关程度，通过计算融合前后影像间的灰度相关系数，可以判断影像融合后的光谱保持程度，若相关系数趋向于1，则说明光谱保持性好，反之则差；SAM是把图像中的每个像元的光谱视为一个高维向量，通过计算两向量间的夹角来度量光谱间的相似性，夹角越小，两光谱越相似，SAM理想值为0[23]；ERGAS主要评价在光谱范围内所有融合波段的光谱质量，考虑光谱变化的整体情况.它的值越接近0，表明在光谱范围内，融合图像的光谱质量越好[24].

a:Landsat8 OLI影像(30 m);b:Sentinel-2影像(10 m);c:蓝波段融合结果(10 m);d:绿波段融合结果(10 m);e:红波段融合结果(10 m);f:近红外波段融合结果(10 m).图2 SFIM融合结果Fig.2 Results of SFIM-based fusion

表1列出了在不同算法下各波段的融合影像质量评价结果，从表中可以看出，ATWT算法的融合效果总体优于SFIM和HPF算法，特别是以Sentinel-2蓝、绿、红3个高分辨率波段进行融合所获得影像的相关系数最高，且SAM和ERGAS值最小，但以近红外波段进行数据融合后，3种算法获得融合影像各波段的相关系数较低，SAM值和ERGAS值高于其他波段融合结果，总体符合目视效果.同时，相比于蓝、绿、近红外波段，在利用Sentinel-2的红波段作为融合时高分辨率参考影像时，融合影像的蓝、绿、红3个波段的相关系数最高，分别达到0.850 405、0.851 331、0.838 229，SAM最低，为5.644 6.

a:Landsat8 OLI影像(30 m);b:Sentinel-2影像(10 m);c:蓝波段融合结果(10 m);d:绿波段融合结果(10 m);e:红波段融合结果(10 m);f:近红外波段融合结果(10 m).图3 HPF融合结果Fig.3 Results of HPF-based fusion

3.2 土地利用分类结果对比

图5 Landsat8影像分类结果Fig.5 Classification result based on Landsat 8 image

结合研究区主要地类分布，分为人工表面、水体、耕地、林地和裸地等5个类别，影像分类结果如图5、图6和图7所示.通过总体分类精度、Kappa系数以及各地类的制图精度、用户精度，对分类结果进行定量分析，结果如表2所示.考虑到分类所用影像空间分辨率的不同，总体而言，3种分类结果的总体精度和Kappa系数之间的关系为：Sentinel-2影像>ATWT融合影像>Landsat影像.ATWT融合影像的总体分类精度及Kappa系数为88.765 7%和0.869 0，低于Sentinel-2的89.693 9%和0.886 0，差距很小，但比Landsat8影像的总体分类精度和Kappa系数提高了6.101 4%和0.056 7.结果表明“红波段+ATWT”融合方法能较好的保留原始影像的光谱信息，虽然在土地利用分类中只能获得接近原始影像分类精度的结果，但该融合方法能为解决Sentinel-2影像中云污染所导致的影像缺失问题提供理论依据，从而为后续融合Sentine-2和Landsat8 OLI数据以获得完整的10 m高分辨率时间序列影像提供了一种有效的方法.

4 结论

本研究结果表明：(1)利用ATWT算法并以Sentinel-2影像的红波段为高分辨率融合参考影像时，融合效果最好.所获得融合影像的B2、B3、B4波段的相关系数最高，分别为0.850 405、0.851 331、0.838 229，同时SAM值最低，为5.644 6.

(2)采用随机森林对整个研究区进行土地利用分类后，融合影像的总体分类精度为88.765 7%，Kappa系数为0.869 0，略低于Sentinel-2影像分类精度，其总体分类精度为89.693 9%，Kappa系数为0.886 0，考虑到空间分辨率的差异，Landsat8影像的总体分类精度为82.664 3%，Kappa系数为0.812 3，远低于之前两种数据的精度.结果表明“红波段+ATWT”融合影像在土地利用分类中可以获得接近原始影像分类精度的结果，能较好的保留原始影像的光谱信息并应用于较大区域的土地利用精细提取.

图6 Sentinel-2原始影像分类结果
Fig.6 Classification result based on Sentinel-2 image

图7 ATWT融合影像分类结果
Fig.7 Classification result based on ATWT fusion image

表2 不同影像分类精度评价表Table 2 Accuracy evaluation of images from different sources

(3)图像融合时选取的Sentinel-2数据与Landsat8数据日期一致，而实际应用中，往往难以获取同一时间的影像，故不同时间Sentinel-2影像与Landsat影像融合在土地利用研究中还有待进一步探讨，以期能得到连续高分辨率时序影像，更好应用于地面种植作物复杂区域的遥感制图.