利用Landsat-8 OLI反演大气气溶胶的可见光谱段地表反射率估算

利用Landsat-8 OLI反演大气气溶胶的可见光谱段地表反射率估算

吕春光1，2，田庆久1，2，王磊1，2，黄彦1，2，耿君1，2

(1.南京大学 国际地球系统科学研究所，南京 210023；2.江苏省地理信息技术重点实验室，南京 210023)

摘要：准确估算地表反射率的贡献一直是遥感反演大气气溶胶光学厚度过程中的重点和难点。为了促进Landsat-8 OLI传感器在地表参数定量化特别是大气遥感领域的应用，本文提出一种利用OLI 1.6μm、2.2μm短波红外谱段数据估算遥感影像可见光地表反射率的方法。该方法依托于MOD04产品地表反射率估算模式，通过光谱归一化和构建新的短波红外植被指数等过程，建立OLI地表反射率估算模式，通过误差分析发现该模式能够有效地降低由于传感器光谱响应不同对估算结果的影响，对应用在OLI遥感影像的计算结果与同时间同区域MOD04产品地表反射率进行比较，表明其结果有较高的相关性和可靠性。

关键词：地表反射率；气溶胶；Landsat-8；光谱归一化；植被指数

doi:10.3969/j.issn.1000-3177.2015.01.008

中图分类号：TP79文献标识码：A

Estimation of Visible Surface Reflectance for Retrieving Aerosol

Optical Depth Using Landsat-8 OLI Data

LV Chun-guang1，2，TIAN Qing-jiu1，2，WANG Lei1，2，HUANG Yan1，2，GENG Jun1，2

(1.InternationalInstituteforEarthSystemScience，NanjingUniversity，Nanjing210023；

2.JiangsuProvincialKeyLaboratoryofGeographicInformationScienceand

Technology，NanjingUniversity，Nanjing210023)

Abstract：Accurate estimation to the contribution of the surface reflectance has been the focus and the difficult point in aerosol optical depth (AOD) retrieval from remote sensing data.In order to promote the applications of Landsat-8 OLI data in quantification of surface physical parameters particularly in the field of atmospheric remote sensing，a method to estimate the optical surface reflectance from OLI images was put forward using OLI 1.6μm and OLI 2.2μm short infrared bands.The method was proposed based on the surface reflectance estimation mode of MOD04 product，by means of processes such as spectral normalization，building a new short infrared vegetation index and so on.The error analysis shows the OLI estimation mode can effectively reduce the error effect caused by spectral response differences.The comparison of the surface reflectance estimation results between the OLI and MOD04 product in the same time and same area represents that the estimation results obtained by OLI estimation mode have high correlation and reliability.

Key words：surface reflectance；aerosol；Landsat-8；spectral normalization；vegetation index

1引言

大气气溶胶是大气中重要的成分之一，它会影响大气辐射状况，对全球和区域气候变化、大气环境质量具有重要的作用和影响，是当前大气遥感研究的重要领域之一[1]。遥感反演大气气溶胶光学厚度(Aerosol Optical Depth，AOD)具有重要的科学和现实意义[2]。有研究指出0.01的地表反射率估算误差将会导致0.1的AOD反演误差，同时将会显著影响气溶胶尺度等物理参数的估算[3-4]，因此，准确估算地表反射率的贡献是AOD反演中的关键环节。

暗目标法是目前最常用的估算方法，它利用浓密植被区红光和蓝光波段地表反射率比较低，且与短波红外地表反射率存在固定的线性关系，从而去除反演中地表对表观反射率的贡献[5-6]。更多的研究表明，在不同观测条件下，可见光波段与短波红外波段地表反射率的比值在一定的范围内波动[7]。中分辨率成像光谱仪(Moderate Resolution Imaging Spectroradiometer，MODIS)C005算法是目前气溶胶反演广泛使用的地表反射率估算方法，它考虑了上述植被覆盖情况和散射角对地表反射率估算的影响[8]，在植被覆盖好的地区具有较高的精度[6]。

但是，由于MODIS气溶胶产品(MOD 04)仅有10km的空间分辨率，从而限制了它的应用范围[9]。2013年2月发射的Landsat-8卫星上搭载的陆地成像仪(Operational Land Imager，OLI)具有30 m空间分辨率，覆盖可见光到短波红外的谱段设置，因此具备利用类似MODIS C005算法的可能性[10]。本文以暗目标法的MODIS地表反射率估算模式为依托，以实测地表反射率数据为基础，考虑传感器的光谱响应关系和谱段设置，尝试建立针对OLI传感器的可见光地表反射率估算模式，为OLI可见光地表反射率的应用提供支持。

2原理与方法

2.1原理

当遥感影像像元(水体像元除外)在2.1μm处的表观反射率满足0.01<ρm2.1<0.25时，该像元可以作为暗目标，在2.1μm处的地表反射率和0.47μm、0.66μm处的地表反射率之间的关系与散射角和植被的茂密程度有关，MODIS地表反射率标准估算模式如下[6-7]：

(1)

(2)

(3)

yint0.66/2.1=-0.00025Θ+0.033

(4)

(5)

(6)

式中，ρs0.66，ρs0.47，ρs2.1分别表示0.66μm，0.47μm，2.1μm波段的地表反射率；ρm1.24，ρm2.1分别表示1.24，2.1μm波段的表观反射率；Θ表示散射角。

图1　OLI和MODIS的不同波段的光谱响应函数与 不同浓密程度植被的光谱对应关系

OLI传感器分别在蓝光波段(0.48μm)、红光波段(0.66μm)和短波红外(2.1μm)3个波段范围与MODIS传感器存在对应关系，但是除红光波段外，每个波段的响应明显要宽于MODIS(图1)。OLI传感器在1.24μm处并没有相应的波段设置，因而不能借助NDVISWIR来建立短波红外与红光反射率之间的关系，但OLI在1.6μm附近的波段位于大气窗口，可考虑与短波红外(2.12μm)结合建立地表反射率的关系。式(4)和式(5)表明散射角Θ在90°～180°变化将会导致MODIS短波红外与红光波段的斜率产生±0.09的差异，同时对截距产生-0.012～0.0105的差异。地表反射率随散射角的变化关系也会随传感器的光谱响应不同而产生差异。

2.2实验与方法

本研究于2013年3月～4月在安徽江苏等地区采集了198组植被冠层光谱数据。实验采用ASD野外光谱仪，其波长范围为350nm～2500nm，其中350nm～1000nm的平均光谱分辨率为3nm，1000nm～2500nm的平均光谱分辨率为10nm，于北京时间11∶00～14∶00之间进行测量，目标光谱测定前、后均立即进行参考板校正，在每个测点采用5次平均值。为了得到OLI传感器地表反射率，并对比与MODIS之间的差别，利用OLI和MODIS的光谱响应函数对实测光谱进行重采样，分别得到MODIS和OLI传感器短波红外2.1μm、2.2μm、1.62μm、1.24μm(仅MODIS)、红光波段和蓝光波段的地表光谱反射率，光谱采样定理如式(7)所示[11]：

(7)

式中，ρs(sensor)为光谱重采样得到的传感器地表反射率，λ为波长，ρs为实测地表反射率，f为传感器光谱响应函数，λmin和λmax为波段起止波长。

3结果与分析

3.1建立OLI地表反射率估算模式

如图2(a)，通过对短波红外(2.1μm、2.2μm)、红光波段和蓝光波段MODIS与OLI地表反射率之间关系可以发现各波段之间分别存在以下关系：

(8)

式中，ρs(MODIS)2.1、ρs(MODIS)0.66、ρs(MODIS)0.47分别为模拟的MODIS地表反射率，ρs(OLI)2.2、ρs(OLI)0.65、ρs(OLI)0.48分别为模拟的OLI地表反射率。如图2(a)所示，其在这两个传感器上有显著的相关性(R2均高于0.99，如图2(a))。将式(8)第2、3式带入式(2)可得到OLI红光波段和蓝光波段的关系：

(9)

图2

(10)

(11)

(12)

利用OLIslopeNDVI′SWIR0.65/2.2替换对应MODISslopeNDVISWIR0.66/2.12，对与每组得到的ρs(MODIS)0.66光谱归一化后的ρs′(OIL)0.65进行比较，使得偏差最小，如式(13)所示：

(13)

由式(13)确定系数后，整理得到以下关系式：

(14)

yint0.65/2.2=-0.000284Θ+0.0352

(15)

将上式(14)、(15)带入式(1)中，并与式(9)、(10)、(11)联立即可得到OLI传感器红光和蓝光波段地表反射率估算模式，Θ为散射角。

针对实测数据光谱重采样得到的OLI 2.2μm反射率，利用上述推导得到的OLI地表反射率模式和直接使用MODIS地表反射率模式对红光和蓝光地表反射进行估算，并分别将估算得到的OLI地表反射率结果与MODIS标准估算模式得到的光谱归一化后的地表反射率结果进行比较，为了表现大小差异，利用下列误差公式进行评价：

(16)

式中，r为地表反射率估算模式误差，ρ′可为利用OLI短波红外波段、OLI估算模式或直接利用MODIS估算模式得到的反射率，ρ为利用MODIS短波红外波段和估算模式计算，并进行OLI光谱归一化后的反射率。图3中，曲线(a)、(b)为对OLI短波红外波段直接利用MODIS标准估算模式估算的误差曲线，其平均误差分别为22%和8.3%；而曲线(c)、(d)为利用建立的OLI估算模式的误差曲线，平均相对误差分别为小于1.5%和1%。这表明考虑传感器光谱响应差异建立OLI估算模式可有效减小MODIS标准估算模式带来的模式误差，采用建立的OLI估算模式估算的地表反射率略低于标准值，而利用MODIS标准模式将会显著高估浓密植被区的地表反射率。曲线(a)、(b)在短波红外小于0.12的范围内，其模式误差r显著大于曲线(c)、(d)，由于浓密植被的2.1μm短波红外反射率主要集中在0.12以下，因此如不考虑传感器光谱响应的差异，直接使用MODIS标准模式的反射率估算关系，其模式本身将会对估计值带来较大误差。

图3　对OLI短波红外数据使用不同地表反射率估算模式 产生的模式误差随2.1μm反射率的变化比较 ((a)直接利用MODIS标准估算模式估算红光波段； (b)直接利用MODIS标准估算模式估算蓝光波段； (c)利用建立的OLI估算模式估算红光波段； (d)利用建立的OLI估算模式估算蓝光波段)

3.2OLI影像地表反射率估算应用与分析

本文选取5景2013年9月1日中国东部地区的Landsat-8 OLI影像数据来进行红光和蓝光波段地表反射率估算。在这5景影像的获取时间范围内天气晴好，各景影像的获取时间、星下点位置和图幅范围如表1所示。选取对应上述OLI影像空间范围的Terra MODIS数据，其时间范围与OLI影像相差在40分钟以内，表1中各序号对应影像的几何位置关系如图4所示。

首先根据影像头文件信息逐像元计算太阳天顶角、方位角、观测天顶角、方位角和散射角，并求得短波红外植被指数；其次，根据植被指数和影像2.1μm波段设定阈值掩膜水体、城镇、云等非暗像元信息，并对2.1μm波段大气校正得到其地表反射率；最后，利用本文建立的OLI估算模式求得红光和蓝光波段地表反射率，结果如图5(a)、(b)所示。

MOD04数据是全球范围使用广泛的MODIS气溶胶产品，它对地表反射率的估算采用了MODIS标准反射率模式[12]。为了对结果进行比较，本文对OLI估算结果进行像元空间重采样，与2013年9月1日MOD04产品的像元进行对应，并从中筛选出非暗像元比例小于10%的样点作为纯像元进行比较。样点OLI和MODIS红光波段地表反射率均值分别为0.0437和0.0485，蓝光波段地表反射率均值分别为0.0279和0.0285。根据像元的经纬度位置进行编号，其中MOD04与OLI地表反射率对应样点的情况如图5(c)所示。从图中可以发现，OLI估算结果在红光和蓝光波段均低于对应点位的MOD04结果，对应曲线的波动趋势具有明显的一致性。对应样点OLI 和MOD04产品短波红外(2.2μm、2.1μm)、红光波段(0.6μm)和蓝光波段(0.4μm)各自的相关关系如图5(d)所示，相关系数R2分别为0.57、0.563和0.538。从图5(d)中散点的分布情况看，MOD04数据某一值通常对应多个分布在周围的不同OLI反射率值，这表明OLI数据对地表反射率的差异更敏感。

表1　所选卫星影像信息

图4　2013年9月1日Terra MODIS 1-2-1波段合成与 Landsat-8 OLI假彩色合成卫星影像位置关系

4讨论

图5(c)表现了在同类型的反射率估算模式下，OLI红光和蓝光地表反射率整体比MODIS低10%和2.3%，导致这种差异的因素很复杂，其中包括传感器的光谱响应、空间尺度效应、地物的反射特性和传感器的观测角度等。由于OLI的观测天顶角很小(最大幅宽小于7°)，本文对OLI地表反射率，利用式(8)进行光谱归一化，发现处理后的红光波段和蓝光波段地表反射率分别提高了3.5%和0.01%，但仍然比MODIS地表反射率低6.6%和2.1%。

图5　OLI影像地表反射率估算结果及其与 MOD04数据的比较

理论上，OLI植被指数的相对增高，会使式(1)短波红外波段2.2μm和红光波段0.65μm反射率之间的斜率相对升高从而导致OLI可见光地表反射率高于MODIS，但实际过程中，散射角的相对大小对斜率有调节作用，并对截距的影响更加显著，从而导致MODIS地表反射率的相对抬升。OLI中所选样点的平均散射角为147.3°，而MODIS对应样点的平均散射角仅为93.7°，根据斜率和截距公式得到各自斜率分别为0.539和0.502，截距分别为-0.005和0.002，而2.1μm平均反射率为0.09，这使得OLI与MODIS相比，受散射角影响的截距抵消了由植被指数和散射角带来的斜率的升高，并进一步将反射率整体降低了0.0037。由此可见，散射角的差异对OLI与MODIS的反射率差异中起决定作用。

以上从地表反射率的影响因素和等式响应机制上对估算结果的差异进行了讨论。此前已有相关的研究表明在传感器对植被观测方向与太阳同向时，其红光反射率有可能较垂直观测时升高[14，18]，这说明MODIS在较大观测角度观测同一目标时，其红光地表反射率有可能比OLI要高。此外，几何定位精度也可能导致两者地表反射率出现差异，MODIS MOD04数据产品和Landsat OLI影像数据本身具备较高的几何校正精度，对二者地表反射率相关性影响较小，由于篇幅和内容所限本文不再对其进行过多的论述。

随着研究的深入，利用红外对可见光地表反射率的估算模式从简单的比值模型发展为引入观测角和植被指数等具有物理意义的复杂模型，也出现了更多不同思路和针对不同传感器估算可见光地表反射率的模式[7，19-23]。不同传感器来探测地表反射率不但受到观测几何、目标特性、波谱响应、瞬时视场和辐射灵敏度的影响，还与传感器的在轨时间、老化程度和卫星的运行状态有关，研究和发展针对不同传感器，并引入更多物理参量提高结果精度和抗干扰因素能力在未来仍具有重要意义和实用价值。

5结束语

本文针对Landsat-8 OLI数据的谱段设置，在MODIS暗目标法估算可见光谱段地表反射率模式的基础上，通过光谱归一化以及构建新的短波红外植被指数等过程，提出了一种新的用于Landsat-8 OLI 2.2μm、1.6μm短波红外谱段数据估算0.65μm红光和0.48μm蓝光谱段地表反射率的方法，并利用该方法与同一时期同类型估算模式的MOD04产品地表反射率进行了比较。本文主要结论有：

(2)利用实测光谱数据的误差分析表明，对OLI谱段直接使用MODIS标准估算模式进行地表反射率计算会产生较大的估算误差，可使红光和蓝光估算结果平均偏高22%和8.3%，而使用建立的OLI估算模式所得红光和蓝光估算误差则分别小于1.5%和1%。这说明考虑不同传感器的光谱响应差异并对估算模式进行改进可有效减小估算模式本身所导致的误差。

(3)对同一地区时间接近的OLI影像地表反射率的估算结果与MODIS产品进行比较发现，结果之间存在较高的线性相关性，而传感器观测角度的不同使得OLI地表反射率略低于MODIS，同时两个传感器地表反射率的差异也是由传感器光谱响应差异、空间尺度效应、观测几何差异、地表植被结构以及卫星运行状态等诸多因素共同作用的结果。

致谢：本文感谢USGS earth explorer网站提供Landsat-8影像数据以及NASA LAADS提供的MODIS数据产品，并对各位审稿专家表示感谢。

参考文献：

[1]DUBOVIK O，SMIRNOV A，HOLBEN B N，et al.Accuracy assessments of aerosol optical properties retrieved from aerosol robotic network (AONET) sun and sky radiance measurements[J].Journal of Geophysical Research-Atmospheres，2000，105(D8)：9791-9806.

[2]TORRE O，BHARTIA P K，HERMAN J R，et al.A long-term record of aerosol optical depth from TOMS observations and comparison to AERONET measurements[J].Journal of the Atmospheric Sciences，2002，59(3)：398-413.

[3]KAUFMAN Y J，TANRE D，GORDON H R，et al.Passive remote sensing of tropospheric aerosol and atmospheric correction for the aerosol effect[J].Journal of Geophysical Research-Atmospheres，1997，102(D14)：16815-16830.

[4]李莘莘，陈良富，陶金花，等.城市与冬季北方亮目标地区气溶胶光学厚度反演[J].中国科学：地球科学，2012，(8)：1253-1263.

[5]KAUFMAN Y J，TANRE D，REMER L A，et al.Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer[J].Journal of Geophysical Research-Atmospheres，1997，102(D14)：17051-17067.

[6]周春艳，柳钦火，唐勇，等.MODIS气溶胶C004、C005产品的对比分析及其在中国北方地区的适用性评价[J].遥感学报，2009，(5)：854-872.

[7]LEVY R C，REMER L A，MATTOO，et al.Second-generation operational algorithm：Retrieval of aerosol properties over land from inversion of moderate resolution imaging spectroradiometer spectral reflectance[J].Journal of Geophysical Research-Atmospheres，2007，112(D13)：1-21.

[8]REMER L A，TANRE D，KAUFMAN Y J，et al.Algorithm for remote sensing of tropospheric aerosol from MODIS：Collection 005[J].National Aeronautics and Space Administration，2006：1-87.

[9]王玲，田庆久，李姗姗.利用MODIS资料反演杭州市500米分辨率气溶胶光学厚度[J].遥感信息，2010，25 (3)：50-54.

[10]MARKHAM B L，STOREY J C，IRONS J R.Landsat data continuity mission-now Landsat-8：Six months on orbit[M].Bellingham：SPIE-Int SOC Optical Engineering，2013.

[11]梁顺林，范闻捷.定量遥感[M].北京：科学出版社，2009.

[12]REMER L A，KAUFMAN Y J，TANRE D，et al.The MODIS aerosol algorithm，products，and validation[J].Journal of the Atmospheric Sciences，2005，62(4)：947-973.

[13]覃文汉，项月琴.植被结构及太阳/观测角度对NDVI的影响[J].环境遥感，1996，(4)：285-290.

[14]张雪红，田庆久，沈润平.冬小麦冠层光谱的方向性特征分析[J].光谱学与光谱分析，2010，(6)：1600-1605.

[15]ZHANG X，YAN G，LI Q，et al.Evaluating the fraction of vegetation cover based on NDVI spatial scale correction model[J].Intnational Journal of Remote Sensing，2006，27(23-24)：5359-5372.

[16]栾海军，田庆久，余涛，等.基于分形理论的NDVI连续空间尺度转换模型研究[J].光谱学与光谱分析，2013，(7)：1857-1862.

[17]ZHANG X H.Study on spatial heterogeneity and scale effect of eucalyptus forest based on high resolution remote sensing[D].Nanjing：Nanjing University，2012.

[18]许民，宜树华，叶柏生，等.植被盖度及太阳/观测角度对疏勒河上游NDVI和SAVI值的影响[J].干旱区资源与环境，2012，(5)：101-107.

[19]TANR E D，DESCHAMPS P Y，DEVAUX C，et al.Estimation of saharan aerosol optical thickness from blurring effects in thematic mapper data[J].Journal of Geophysical Research：Atmospheres (1984-2012)，1988，93 (D12)：15955-15964.

[20]盛莉，黄敬峰，郭瑞芳，等.综合利用环境星CCD和红外数据反演大气气溶胶光学厚度[J].红外与毫米波学报，2013，(3)：277-282.

[21]王中挺，厉青，王桥，等.利用深蓝算法从HJ-1数据反演陆地气溶胶[J].遥感学报，2012，(3)：596-610.

[22]HAU N C，TSAY S C，KING M D，et al.Deep blue retrievals of asian aerosol properties during ACE-Asia[J].IEEE Transactions on Geoscience and Remote Sensing，2006，44(111)：3180-3195.

[23]KAUFMAN Y J，WALD A E，REMER L A，et al.The MODIS 2.1-mu m channel-correlation with visible reflectance for use in remote sensing of aerosol[J].IEEE Transactions on Geoscience and Remote Sensing，1997，35(5)：1286-1298.