Sentinel-2A与Landsat 8 OLI逐像元辐射归一化方法研究

徐玉雯，张 浩，陈正超，景海涛

(1.河南理工大学测绘与国土信息工程学院，河南 焦作 454003；2.中国科学院遥感与数字地球研究所数字地球重点实验室，北京 100094)

随着在轨运行卫星数量和种类的增加，越来越多的遥感应用倾向于综合应用多源遥感数据开展高频次的应用监测，这就要求不同传感器数据解译具有一致性，通过辐射归一化使数据在空间和时间上保持连续性，使定量研究的重点聚焦在地表的变化而不过分考虑传感器的差异。目前常用的辐射归一化方法可分为绝对辐射归一化(absolute radiometric normalization,ARN)及相对辐射归一化(relative radiometric normalization,RRN)[1]。绝对辐射归一化是将图像的DN值转换成地表反射率。相对辐射归一化是以同区域多时相相似波段灰度值之间满足线性关系为前提，主要分为基于分布的相对辐射归一化方法和基于像元对的辐射归一化[2]。基于分布的相对辐射归一化方法简单、客观，考虑全局分布特征，归一化结果受主观因素影响较少，在多源遥感器之间的辐射归一化更显优势[3]。基于像元对相对辐射归一化方法主要原理是从待归一化的影像对中选取部分或全部的像元对作为样本点回归得到校正系数，而样本点质量的好坏直接影响归一化校正的结果[4-9]。上述相对辐射校正方法，不同影像相似波段之间的变换均属于线性变换，即同一波段所有像元具有统一的变换系数，不同的是选取伪不变特征点的方式及回归方程的拟合方法有所差别，并且这些方法多应用于同类型传感器多时相影像之间的归一化，没有考虑不同传感器对相同地物的成像差异及不同地物之间的光谱差异。文献[10]考虑地物差异利用表观反射率光谱库，通过回归分析与线性组合变换实现了SPOT 5与THEOS之间的辐射归一化，但是该归一化是基于表观辐射亮度进行的，没有在地表反射率层面上实现辐射归一化[10]。本文考虑不同传感器光谱响应差异对地物反射率的影响，以两景几乎同时成像的Sentinel-2A和Landsat 8 OLI数据为例，基于光谱库构建两传感器VNIR相似波段之间的匹配转换模型，以实现逐像元归一化。

1 研究区数据及预处理

本文采用准同步过境Sentinel-2A及Landsat 8 OLI为数据源，以内蒙古中东部达里诺尔湖附近地区为研究区，两景影像在2017年7月17日成像时间分别为3:05:41、2:58:44，相差7 min。Sentinel-2是全球环境与安全监测计划(GMES)下的高空间分辨率卫星，包括A、B两星影像，单景成像幅宽与用户分发幅宽分别为290、109.8 km，光谱范围为0.4～2.4 μm，包含13个波段，空间分辨率分为10、20、60 m 3个等级[11]，其Level-1C数据是经正射校正和亚像元级几何精校正后的表观反射率产品。Landsat 8 OLI陆地成像仪成像宽幅为185 km，包括9个波段，空间分辨率为30 m，其中全色波段15 m，L1T产品是使用地面控制点和数字高程模型数据进行精确校正后的DN值数据[12]。同作为提供全球免费观测服务的新型中高分辨率传感器，两者在定量遥感方面应用广泛，波段设置也较为相似，VNIR波段光谱响应函数如图1所示，波段参数对比见表1。

表1 Sentitle-2A与Landsat 8 OLI可见—近红外波段对比

数据预处理主要包括大气校正、无效值掩模、空间重采样、裁剪等。Sentinel-2A影像采用欧空局官方提供的L2A_AtmCorr(Sentinel-2 atmospheric correction)算法反演地表反射率，其中气溶胶和水汽的反演分别基于改进的DDV(dense dark vegetation)算法及APDA(atmospheric pre-corrected differential absorption)算法，校正得到的地表反射率产品以整型存储，放大系数为10 000，有效值范围0～10 000[13]。Landsat 8 OLI影像采用NASA官方提供大气校正算法进行大气校正，该算法是基于MODIS影像获取相应大气参数经反演得到地表反射率[14]，放大系数为10 000，地表反射率可能值为-2000～12 000，有效值范围0～10 000。无效值对影像归一化及后续分析并无意义，需要对Sentinel-2A与Landsat 8图像上无效反射率值进行掩模处理。然后将Sentinel-2A的蓝、绿、红波段10 m与近红外20 m的空间分辨率采用像素聚合方法重采样到与Landsat 8相同的空间分辨率，根据二者同时刻成像重叠情况，从配准好的两影像中裁剪出600×600像元大小的图像作为研究区，如图2所示。研究区中包含了浓密植被、草地、高亮白沙地、水体等主要地物。

2 逐像元归一化方法

2.1 光谱匹配因子库建立

光谱库是由光谱成像仪在实验室或野外对不同地物按照严格规范的流程取的各类地物反射光谱数据的集合，具有多样性、稳定性的特点，对遥感图像的解译、影像地物的匹配识别具有重要的作用。ENVI中内嵌的光谱库包含了国际上常用的5种光谱库：USGS、JPL、JHU、IGCP-264、ASTER，近2000条光谱，但并不是所有的地物光谱都涵盖VNIR波长范围的反射率值。对经过筛选的反射率光谱整合、重采样，通过与Sentinel-2A和Landsat 8 OLI的VNIR波段的光谱响应函数按照式(1)进行等效计算，形成不同地物类别的等效光谱库。同一地物因传感器光谱特性差异反映出的波段等效值也不同，两传感器等效光谱库对应地物之间的比值则组成了光谱匹配因子库。

(1)

式中，ρj为某一地物波段j的等效反射率；λ1、λ2为该波段的起止波长；Rj(λ)为波段j的响应函数，ρ对应该波段光谱库实际反射率值，Rj(λ)和ρ是间隔1 nm的不连续的值，对其积分通常用累加值替代。

另外，由于ENVI自带光谱库中缺少清澈水体的标准光谱，试验中又无实测类似的光谱数据。又因达里诺尔湖地处内蒙古中东地区，水体受污染和破坏的程度较小，查询往年同地区相同季节的Hyperion高光谱数据，下载得到2012年7月27日一景影像。经过坏波段剔除、坏线条纹去除、减弱smlie效应、大气校正等处理后获取可见近红外的48个有效波段，在达里诺尔湖上抽样100个水体的像元点平均以获取连续的水体光谱，录入光谱库，与光谱库其他光谱同作处理。

2.2 逐像元归一化模型

传感器成像过程中，每个像元对应的实际地物往往是由一种或多种地物混合而成的，反映出的地物光谱反射率值也是一种或多种地物光谱的组合。通过完全约束最小二乘约束对Sentinel-2A图像进行线性分解获取典型地物的丰度图像，丰度图像客观地表示了每种地物在每个像元中所占的比例，而每个混合像元的光谱反映值则是端元光谱与丰度的线性组合[14-17]。

最小二乘约束的原则是像元内每种端元的比例值范围在0～1之间，并且和为1，即

∑wi=1wi∈[0,1]

(2)

基于光谱解混原理，利用光谱库中匹配到的标准光谱去表示影像像元地物的光谱，通过统计分析和多元线性回归分析，实现影像实际光谱到光谱库光谱的转换。假设ρn{ρn1,ρn2,…,ρni}为归一化时所用目标影像Sentinel-2A各波段各像元点的光谱集合，sn{s1n,s2n,…,skn}表示从光谱库中经匹配分析获取的像元代表光谱在波段n处的等效值，设目标影像像元多元线性回归表达式为

ρni=a0i+w1is1ni+w2is2ni+…+wkiskni

(3)

式中，i为像元点数；n为波段数，参考影像与目标影像的相似波段一一对应；wki为权重，表示光谱库中某一匹配光谱占目标影像实际像元光谱的比重；k表示典型光谱的数量。各代表光谱经匹配因子gkn调整之后，按像元顺序组合，则得到归一化之后各波段各像元点的光谱集合ρN{ρN1,ρN2,…,ρNi}，即

lkNi=gknskni

(4)

ρNi=b0i+w1il1Ni+w2nl2Ni+…+wkilkNi

(5)

3 数据结果处理与分析

3.1 典型地物解混

根据影像中地物类别，结合几种经典的选取端元的方法及目视判断，选取植被、草地、高亮白沙地、水体4种典型地物作为合适的待选端元，综合采用最小距离及光谱角法在等效光谱库中匹配和典型地物光谱最相似的标准地物光谱，并将其等效值作为输入对Sentient-2A影像进行解混，得到的丰度图像如图3所示。从解混结果上看，各典型地物的代表性可以较好地体现出来。

3.2 归一化结果对比

基于4种典型地物的比例图像及依据光谱库计算的匹配因子库，按照式(5)逐像元校正组合之后获取归一化之后的Sentinel-2A影像，同原始影像及采用HLS[17]线性参数归一化后影像，分别与同时刻的Landsat 8参考影像进行散点回归分析，如图4所示。散点图中实线为线性拟合线，虚线为y=x,拟合方程、R2、RMSE统计见表2。

表2 相关性分析统计

从图表中可以看出，虽然Sentinel-2A与Landsat 8 OLI的波段设置相似，选取的影像对成像时间仅相差7 min，大气条件接近，地物类型也相同，但是经过各自官方大气校正算法反演的反射率影像却具有差异。对比原始的Sentinel-2A反射率影像，发现未归一化的Sentinel-2A VNIR波段反射率比Landsat 8 OLI偏高，蓝波段最为明显，且从蓝波段到近红外波段这种偏高的趋势依次递减，相关性逐渐增加，各波段的R2依次为0.805 9、0.935 2、0.961 1、0.991 2。即随着波长的增加，两传感器影像反演得到的反射率更加相近，受大气校正方法差异及辐射响应差异的影响越来越小。但是从影像中各地物类型分开来看，并不是所有地物的各波段的Sentinel-2A反射率值均比Landsat 8 OLI高，比如植被区，在红波段Sentinel-2A的反射率值比Landsat 8 OLI低，而其他3个波段则相反，4个波段在不同类型的地物区呈现不同的高低值。

HLS计划下的光谱通道归一化系数是通过拟合模拟的Landsat 8和Sentinel-2地表反射率值得到的。从中间的散点图可以看出，同原始影像一样，经HLS校正系数归一化的影像R2随波长的增加而增大，相比于原始影像，蓝波段相关性有细微的提高，红、近红外波段保持不变，绿波段降低。其次，就RMSE而言，蓝、红波段的RMSE相对于原始影像各减少了1.1%、2.7%，绿、近红外波段保持不变。综合来看，经HLS校正系数归一化之后，蓝、红波段效果有所改善，但是并不明显，近红外波段无变化，绿波段引入了额外的误差。这可能由HLS系数针对的是相对一致的大气校正方法，而本文采用的是各自官方的大气校正算法所致。

逐像元归一化后的Sentinel-2A VNIR 4个波段从R2、RMSE、拟合系数3个方面看，都更接近于参考影像，辐射一致性明显增强。相比原始影像R2，除红波段稍差之外，其他3个波段均有所提高，蓝波段提高的最为明显，直接从0.805 9增至0.813 2。从RMSE分析，蓝、绿、近红波段RMSE各降低9.7%、1.5%、1.4%，近红外保持不变。从拟合方程上看，经逐像元归一化后各波段的slope更接近于1，蓝、绿、红、近红外波段的slope分别由原来的1.131 1、1.182 6、1.181 5、1.055 2降低到1.040 0、1.178 9、1.149 0、1.054 1。相对于一个波段所有像元均采用同样的归一化校正系数的方法，本文提出的通过匹配因子调整的逐像元归一化方法可以一定程度上减少传感器之间的差异，增强多源传感器之间的辐射一致性，有利于遥感影像的综合利用。

4 结 语

Sentinel-2A与Landsat 8 OLI均是新一代中高分辨率传感器，提供全球观测服务，波段设置相似，数据也可免费获取，且均具有大幅宽、高精度、高空间分辨率的特点，在定量遥感方面应用广泛。为了方便用户更加高效合理地利用影像数据，两者在众多大气校正算法的基础上发展了各自的反射率反演算法，但即便是同时刻成像的影像，得到的反射率数据差异也比较大。因此，对其综合应用时，需要考虑两者大气校正算法的特点及光谱响应的差异，将其归一化至同一辐射基准。

基于光谱库的逐像元归一化方法可以有效地消除部分光谱响应带来的差异，相对于传统简单的相对辐射归一化所有波段均采用一个套转换系数，本文方法考虑地物的类型及传感器波段设置特点，经归一化后的影像与目标影像各波段之间的相关性有所提高，RMSE降低，辐射一致性增强。更重要的是，与传统方法相比，该方法为多源中高分辨率遥感图像高精度辐射归一化提供了新思路。

然而，在试验过程中也发现，当地物类型复杂时，典型地物的选择、组合、光谱匹配精度等因素在一定程度影响归一化效果。今后研究工作考虑在传统辐射归一化的基础之上再应用本文方法，结合两者的优势，在其他多源中高分辨率反射率影像上验证。