0 引言

近年来, 卫星遥感技术发展迅速, 众多新型卫星相继发射, 遥感对地观测能力逐步提高.然而, 由于受到云、雨天气的影响及回归周期的限制, 单颗光学遥感卫星在进行长时间序列和高质量的无云影像获取中皆显得力不从心, 为此, 实现多源遥感数据的互补应用则显得尤为重要.徐涵秋等^[1]利用归一化植被指数(NDVI)将GF-1 PMS与ZY-3 MUX传感器在植被观测能力上进行了对比, 发现ZY-3 MUX的NDVI信息量和信号总量强于GF-1 PMS. Van Leeuwen等^[2]对MODIS和AVHRR影像的NDVI进行对比, 发现两者有较好的线性相关性, 通过线性方程可相互转换, 但NDVI指标仅涉及到近红外和红波段, 没有对其它波段信息进行对比.贾玉秋等^[3]利用GF-1与Landsat-8多光谱遥感影像进行了玉米LAI反演对比, 发现GF-1影像的蓝、绿、红和近红外波段表观反射率(Top of the Atmosphere Reflectance, TOA反射率)与Landsat-8有显著的线性关系, 但所获得的线性关系方程由于选取的样本点较少, 转换精度有待提高.苏涛等^[4]将IRS-P6 LISS-3和Landsat-5 TM两种多光谱传感器影像进行了数据对比, 也针对对应波段的TOA反射率建立了转换方程, 但该研究并没有考虑各波段的地表真实反射率(Bottom of the Atmosphere Reflectance, BOA反射率).当前在轨运行的代表性资源卫星中, Landsat-8是美国国家航空航天局(NASA)于2013年2月11日发射的第八代陆地资源卫星^[5]; 2015年6月23日发射的哨兵2A号(Sentinel-2A)是欧空局(ESA)"哥白尼"计划下多光谱成像任务中的首发星^[6]; 高分一号(GF-1)卫星是我国高分辨率对地观测系统重大专项天基系统中的首发星, 于2013年4月26日发射^[7]. 3颗卫星凭借数据精度高、更新速度快、成像波段多等优势, 广泛应用于精细农业、污染监测、土地资源管理、公共安全等领域.可以预见, 这3种卫星数据在当前以及今后相当长时期内, 将成为最具广泛应用前景的地表资源环境遥感监测数据源.目前3颗卫星运行状况良好, 数据质量高且可免费获取, 是今后实现卫星长时间序列数据集构建重要的互补性数据源.因此, 开展上述卫星数据间的光谱信息对比及定量分析转换研究具有重要意义.基于此, 本研究将GF-1与当前国际上使用广泛的Landsat-8以及Sentinel-2A进行对比, 定量分析GF-1与另外两颗卫星在蓝、绿、红、近红外波段TOA反射率以及BOA反射率中的差异和相关性, 以期为提升我国GF-1数据的应用价值、拓展数据的应用领域、增强多源多光谱遥感数据在定量遥感研究中的互补提供技术支持.

1 数据

为了准确量化不同卫星多光谱传感器各波段之间的信息差异, 一个重要的前提条件是选取同日过境、太阳高度角和方位角接近且有重叠覆盖范围的无云卫星影像.这样可以保证传感器在获取遥感数据时大气影响的差异最小, 并且获得的地表反射能量最为接近.由于不同卫星的重访周期和幅宽不同, 因此难以获取一个地区内的同一天、同时过境的三颗卫星影像.本研究获取了2014年9月8日同天过境研究区1(123$^{\circ}$02$'\sim$123$^{\circ}$44$'$E, 41$^{\circ}$11$'\sim$41$^{\circ}$37$'$N)的GF-1 WFV2及Landsat-8 OLI遥感影像, 以及2016年8月22日覆盖研究区2(122$^{\circ}$48$'\sim$123$^{\circ}$47$'$E, 41$^{\circ}$10$'\sim$41$^{\circ}$35$'$N)的GF-1 WFV3及Sentinel-2A MSI影像.研究区主要土地覆盖类型包括耕地、林地、水体、裸土和不透水面(建筑物、道路等), 具体地理范围如图 1所示.

GF-1卫星搭载了一台2 m分辨率的全色相机、一台8 m分辨率的多光谱相机(PMS)和四台16 m分辨率多光谱相机(WFV).基于卫星观测模式和对比影像的空间分辨率信息特征, 共获取两景GF-1多光谱16 m影像.研究区1选取的GF-1 WFV2影像成像时间为2014年9月8日, 研究区2选取的GF-1 WFV3影像成像时间为2016年8月22日. Landsat-8卫星携带有OLI和TIRS两个传感器: ①陆地成像仪OLI, 包括9个波段, 空间分辨率为30 m(全色波段为15 m), 成像幅宽为185 km; ②热红外传感器TIRS, 包括2个热红外波段, 空间分辨率为100 m.本研究获取的Landsat-8 OLI影像包括L1T级产品和在其基础上经过大气校正的L2级产品, 成像时间为2014年9月8日. Sentinel-2A卫星携带一个多光谱传感器MSI, MSI包括13个波段, 空间分辨率有10 m、20 m、60 m 3种, 成像幅宽290 km.本研究获取的Sentinel-2A MSI影像成像时间为2016年8月22日.上述3颗卫星的具体参数及要对比的对应波段信息如表 1及图 2所示.同时, 课题组于2016年9月15-17日对研究区进行了为期3 d的实地调查.通过田间观测和走访咨询, 获得了大量研究区的地物类型(耕地、林地、水体、城镇等)、分布规律及植被的物候特征等信息(见图 2), 这为提高遥感影像的地类识别, 以及对比样本点选择的精度提供了有效保证.

2 方法 2.1 数据处理

对于GF-1的WFV2和WFV3影像, 处理步骤包括:正射校正, 几何校正, 辐射定标和大气校正.首先, 参考空间分辨率为30 m的DEM数据对GF-1影像进行正射校正, 减小因地形起伏对数据造成的影响; 然后, 分别以经过严格校正的Landsat-8以及Sentinel-2A影像作为参考^[8], 对同时期GF-1数据进行几何校正; 接下来对校正后的数据进行辐射定标, 定标参数从中国资源卫星应用中心获取(http://www.cresda.com/CN/Downloads/dbcs/10506.shtml); 最后, 通过ENVI5.3软件提供的FLAASH方法进行大气校正.具体定标公式如下.

GF-1辐射定标公式:

$ L_{\lambda}=Gain\cdot DN. $

(1)

式中, $L_\lambda$为辐射亮度, 单位为W/(m$^{2}\cdot $sr$\cdot \mu $m); $Gain$为增益; $DN$为亮度值. GF-1辐射亮度转换为TOA反射率计算公式为^[9]

$ \rho_{\lambda}=\dfrac{\pi\cdot L_{\lambda}\cdot d^2}{E_S\cdot \cos(\theta_{SE})}=\dfrac{\pi\cdot L_{\lambda}\cdot d^2}{E_S\cdot \sin(\theta _{SZ})}. $

(2)

式中, $\rho_{\lambda}$为TOA反射率; $L_{\lambda}$为辐射亮度; $d$为日地天文单位距离; $E_S$为大气层外波段太阳辐照度, 单位为W/(m$^{2}\cdot \mu $m); $\theta_{SE}$为太阳天顶角; $\theta_{SZ}$为太阳高度角, 且满足$\theta_{SE}=90^{\circ}-\theta_{SZ}$.

对于Landsat-8 OLI的L1T级产品, 主要通过美国地质调查局(USGS)提供的定标参数(可从影像头文件中获得)对其进行辐射定标处理, 得到大气上行辐射亮度, 再计算TOA反射率. L2级产品为经过大气校正之后的BOA反射率数据, 经USGS数据分发网站(http://earthexplorer.usgs.gov/)可以免费下载.本研究采用的具体计算公式如下.

Landsat-8 OLI大气上层辐照度计算公式^[10]:

$ L_{\lambda}=M_L\cdot Q_{CAL}+A_L. $

(3)

式中, $L_\lambda$为辐射亮度, 单位为W/(m$^{2}\cdot $sr$\cdot \mu $m); $M_L$为辐射亮度增益; $A_L$为辐射亮度偏移; $Q_{CAL}$为影像DN值. Landsat-8 OLI辐亮度转换为反射率计算公式:

$ \rho_\lambda=M_\rho\cdot Q_{CAL}+A_{\rho}/ \cos(\theta_{SZ}). $

(4)

式中, $\rho_{\lambda}$为TOA反射率; $M_{\rho}$为反射率增益; $A_{\rho}$为反射率偏移; $Q_{CAL}$为像元$DN$值; $\theta_{SZ}$为太阳高度角. Sentinel-2A MSI影像是从ESA网站(https://scihub.copernicus.eu/dhus/#/home)下载到的, Sentinel-2A MSI L1C级产品为经过定标之后的TOA反射率, 为方便数据的存储, 数值被扩大了10 000倍^[11].主要数据处理步骤是进行大气校正.本研究采用的是欧空局开发的开源遥感数据处理软件SNAP上用于处理Sentinel-2A L1C级产品的工具箱SEN2COR来计算BOA反射率^[12].

2.2 样本数据选择

为了对比两组传感器在蓝、绿、红、近红外波段TOA反射率以及BOA反射率的差异, 参考同时期高分2号(GF-2)的高清影像(0.8 m)和2016年9月的实地调查数据, 选取了耕地、林地、水体、裸土和不透水面等类型用于反射率定量对比研究.样点数据的选择是本研究中最为关键的步骤, 后期数据的统计分析都在此基础上进行.为此, 样点的选择应满足以下原则.

(1) 选取无云区域, 尽可能减少复杂大气环境对结果造成的影响.

(2) 样点应选择均质区域, 并尽可能选取均质区域中央位置, 以确保选择纯像元.

(3) 应在一定的波长范围内选取样本, 即高低反射率地物都应涉及, 这样才能客观地对比不同传感器在整个波谱范围的性能.

(4) 选取地每类地物样本数目不少于300个, 以确保不同传感器对应波段间数据的可信度.

GF-1 WFV2与Landsat-8 OLI TOA反射率和BOA反射率的对比研究共选取了1 903个纯样点, GF-1 WFV3与Sentinel-2AMSI的对比则选取1 527个纯样本点.对于每组对比图像, 随机选取每类地物纯样本点中2/3用于不同传感器光谱信息对比及回归建模, 其余的1/3用于转换结果的精度验证.

2.3 回归模型建立与检验

本研究采用统计线性回归方法, 对不同传感器各个对应波段的TOA反射率及BOA反射率进行逐一对比.利用决定系数$(R^{2})$和均方根误差(RMSE)这两个评价指标, 分别进行对应波段光谱反射率回归分析.然后利用回归得到的转换方程进行TOA反射率及BOA反射率的模拟, 线性回归方程的显著性检验是借助$F$检验来实现的.

决定系数$R^{2}$均方根误差RMSE计算公式为^[13]

$ R^2=\dfrac{\sum^n_{i=1}(\widehat{y}_l-\overline{y})^2} {\sum^n_{i=1}(y_i-\overline{y})^2}, $

(5)

$ {\rm RMSE}=\sqrt{\dfrac{\sum^n_i(\widehat{y}_l-y_i)^2}{n}}. $

(6)

式中, $n$为验证样点数量; $y_i$为样点实测值; $\widehat{y}_1$为样点估测值; $\overline{y}_l$为样点平均值.

$F$检验计算公式为

$ U=\sum\nolimits^n_{i=1}(\widehat{y}_l-\overline{y})^2, $

(7)

$ Q=\sum\nolimits^n_{i=1}(y_i-\widehat{y_l})^2, $

(8)

$ F=\dfrac{U}{Q/(n-2)}. $

(9)

式中, $Q$为误差平方和或剩余平方和; $U$为回归平方和; $n$为样本数.

3 结果与分析

在选择的样本点处, 分别提取了3颗卫星各个波段的TOA及BOA反射率数值, 得到了两组卫星对应波段的光谱信息对比结果(见图 3).

由图 3可见, 对于TOA和BOA反射率, GF-1 WFV与Landsat-8 OLI和Sentinel-2A MSI数据在蓝(B)、绿(G)、红(R)和近红外(NIR)波段都不同程度地存在一定偏差.而且, 对比卫星对应波段的光谱值的相关系数$r$均大于0.96, 这表明各个对应波段反射率都具有很强的线性相关性.横向对比结果显示:在TOA反射率上, GF-1 WFV2各个波段反射率数值普遍高于Landsat-8 OLI, 拟合直线与1：1线几乎平行(图 3(a1)-(d1)); GF-1 WFV3各个波段反射率数值普遍高于Sentinel-2A, 拟合直线随着反射率数值增加而逐渐偏离1：1线(图 3(a2)-(d2)); 对于BOA反射率, GF-1 WFV2与Landsat-8 OLI各个波段的差异不大, 分布在1：1线两侧(图 3(a3)-(d3)); GF-1 WFV3各个波段BOA反射率数值大部分高于Sentinel-2A MSI, 并随着反射率数值增加差异逐渐增大(图 3(a4)-(d4)).通过纵向对比可发现, 与TOA反射率对比结果相比, 经过大气校正之后的BOA反射率数值的差异降低明显, 一致性明显提高.

为了分析产生上述波谱信息差异的原因, 本研究绘制了3颗卫星对应波段的光谱响应曲线(见图 4).其中, GF-1WFV光谱响应曲线由中国资源卫星应用中心提供的光谱响应函数进行绘制(http://www.cresda.com/CN/Downloads/gpxyhs/index.shtml); Landsat-8 OLI的光谱响应曲线则通过美国地质调查局(USGS)提供的光谱响应函数进行绘制(https://landsat.usgs.gov/using-usgs-spectral-viewer); Sentinel-2A光谱响应值可直接从下载到的影像参数信息中获取并绘制成图.

由图 4可以看出, GF-1 WFV2与Landsat-8 OLI两传感器的蓝与绿波段波谱带宽设置较为接近, 但在红和近红外波段, 波谱带宽差异显著(见图 4).对于GF-1 WFV3与Sentinel-2A MSI而言, 除了近红外波段, GF-1的波谱带宽明显要比Sentinel-2A宽.

波谱范围越宽, 则传感器获得的能量越多.由图 4(a)可知, WFV2各个波段的波谱范围比OLI要宽, 在蓝、绿波段, 两者光谱响应值差异不大, 在红波段中心波长位置WFV2光谱响应值稍低于OLI.在近红外波段, 只有在很窄的范围内WFV2光谱响应值小于OLI, 光谱响应值越大, 传感器接收到的信号能量越强.结合图 3(a1)-(d1)不难发现, WFV2在蓝、绿波段TOA反射率比OLI要大, 在红、近红外波段部分值要小于OLI, 这与光谱响应曲线体现出来的差异一致, 波谱范围和光谱响应值共同影响反射率的大小.由图 4(b)可见, WFV3各个波段的波谱范围大于Sentinel-2AMSI, 在近红外波段波谱范围差异最小.各对应波段WFV3光谱响应值比MSI要高, 在蓝、绿、红波段差异最大, 在近红外波段差异最小.通过图 3(a2)-(d2)容易发现, WFV2各个波段的TOA反射率数值要大于MSI反射率数值, 蓝、绿、红波段反射率的差异最大, 近红外波段反射率差异最小, 两者结果一致.由图 3(a3)-(d3)和图 3(a4)-(d4)发现, BOA反射率数值比TOA反射率数值低, 由于大气校正, 减弱了大气的影响, 降低了反射率.

基于对两组数据的波段信息对比可知, 构建GF-1与其他两颗卫星对应波段的光谱信息转换方程, 便可实现不同卫星之间对应波段的相互转换.通过对GF-1 WFV与Landsat-8 OLI和Sentinel-2A MSI影像各个波段反射率分别进行回归分析, 获得了WFV2与OLI以及WFV3与MSI传感器在蓝、绿、红、近红外波段的反射率转换方程(见表 2).

为了检验转换方程的效果, 本研究将Landsat-8 OLI和Sentinel-2A MSI进行了波段转换, 并将转换结果分别于GF-1 WFV2及GF-1 WFV3进行了对比验证.其中, 选择了总体样本点中余下的1/3作为验证数据用来进行光谱转换效果的验证, 结果如图 5所示.

通过对比由转换方程计算的GF-1 WFV2和WFV3的各个波段TOA反射率与另外两颗卫星数据的对应波段反射率值发现, 两组对比数据的一致性显著提高, 拟合直线与1：1线几乎重叠(见图 5).统计分析显示, 对于GF-1 WFV2与Landsat-8 OLI两传感器在可见光(蓝、绿、红)波段, RMSE都低于0.008, 在近红外波段RMSE < 0.015, $R^{2}>$0.93.通过一元线性回归模型$F$检验得到, 在置信度水平$\alpha=0.05$下, $F\gg F_{0.05}(1, 632)$, 说明回归方程在此水平上是显著的.对于GF-1 WFV3与Sentinel-2A MSI两传感器, 可见光(蓝、绿、红)波段的RMSE < 0.007, 在近红外波段RMSE最大值为0.0181, $R^{2}>$0.92. $F$检验结果显示在置信度水平$\alpha=0.05$下, $F\gg F_{0.05}(1, 509)$, 这表明回归方程在此水平上是显著的, 所建立的光谱转换方程的转换效果较好.

4 结论

本文通过选取同日过境的卫星影像, 对比分析了GF-1 WFV2与Landsat-8 OLI以及GF-1 WFV3与Sentinel-2A MSI对应波段在TOA及BOA反射率上的差异.通过统计回归拟合得到了GF-1与另两颗卫星的波段光谱信息转换方程, 并利用大量的样本点对光谱转换方程进行了验证.研究表明, 多光谱数据GF-1与Landsat-8以及GF-1与Sentinel-2A的各个对应波段在TOA以及BOA反射率上具有很好的线性相关性, 对于少量云层覆盖, 地物类型丰富的地区, 可以对不同传感器数据采用线性光谱转换.在进行光谱转换时, 两景影像时间间隔不宜过长, 否则影响转换结果.通过转换方程计算得到的GF-1转换值与观测值对比验证得出, 本研究获得的GF-1 WFV2与Landsat-8 OLI以及GF-1 WFV3与Sentinel-2A MSI转换方程在TOA和BOA反射率中能够实现数据间的高精度转换, 这将为多源遥感数据在农作物生长状态的长期定量监测、自然灾害持续定量监测等领域提供技术支持.

另外, 本研究显示, 所选择的3个卫星的过境时间存在差异(约0.5$\sim$1 h)(见表 1), 因此如果研究区域有云移动的情况, 可以用无云区域的影像来进行互相弥补.本研究主要针对的是同日过境的卫星影像进行波段信息对比及转换, 而对于相隔一定天数比如间隔1 d、2 d、3 d天等等所获取数据间的光谱关系则需要进一步分析, 这也是今后研究的重点.总之, 时间间隔越大, 误差越大, 原因是由于大气状况以及地表状况将随时间发生明显变化, 必将导致各传感器波段的光谱信息发生变化.

虽然本研究获得了GF-1与Landsat-8以及Sentinel-2A在蓝、绿、红、近红外波段的TOA以及BOA反射率的转换方程, 但由于数据获取限制, 对于分类型的地表反射率差异分析及转换的研究依然缺乏, 尚需进一步检验及深入探讨.由于3颗卫星空间分辨率各不相同, 过境时间也略有偏差, 再加之不同大气校正算法的差异也在一定程度上对结果造成偏差.本研究只对同日过境的卫星影像进行了波段光谱信息对比以及转换, 但现实应用中, 大多数遥感数据的获取时间都存在一定差异, 对于相隔不同天数的数据光谱之间的定量关系有待进一步探讨.因此, 从卫星遥感数据协同应用的角度讲, 这将成为研究的一个出发点, 不仅有助于今后更好地利用已有的、免费共享的多光谱遥感数据, 而且对于加强多源遥感卫星数据在长期、连续的地表参数定量反演的应用潜力都具有重要价值.

致谢: 感谢中国资源卫星应用中心、美国地质调查局以及欧空局为本研究提供遥感数据.