1 引言

光学卫星影像是遥感应用的重要数据源之一。随着国内外航天技术的不断发展,可用的卫星影像种类不断增加,时间分辨率、空间分辨率和光谱分辨率也在不断提高,为遥感数据的定量化应用提供了基础。Sentinel-2影像是新一代多光谱影像,共有13个光谱波段,其中包括3个红边波段,是植被长势监测的理想数据源,在全球环境变化、农作物生长监测、森林监测、灾害应急与救援等方面具有重要意义^[1]。

在太阳-大气-目标-大气-传感器的辐射传输过程中,由于大气分子、气溶胶的散射以及臭氧、水汽等气体的吸收,卫星接收的信号不能准确地反映地表物理特征,从而造成影像的清晰度和对比度下降以及反射率、辐亮度等相关物理量出现偏差,降低了卫星影像的质量和可用性^[2]。大气校正的目的是通过一定的手段消除或减少大气分子、气溶胶的散射和臭氧、水汽等气体吸收对地物反射率的影响^[3],是进行地表参数定量分析的前提。

大气校正的方法主要有经验线性法、物理模型法和介于经验方法和物理模型之间的半经验物理模型法。经验线性法是基于影像统计量或地面实测的实验方法,包括直方图调整法、暗物体扣除法、固定目标法等^[4],其优点是简单、易于操作,缺点是依赖各种假设或实测数据,对影像的要求较高,容易出现过度校正。基于大气辐射传输模型的理论方法即物理模型法,如6S模型、FLAASH模型、MORTRAN模型、ACTOR模型,以数学、物理学为基础,通过研究太阳辐射在大气介质中的传输过程,得到大气介质的吸收和散射特性。相比经验线性法,基于大气辐射传输模型法具有明确的物理意义,精度更高,是目前使用最广泛的大气校正方法。介于经验方法和物理模型之间的半经验物理模型法,如大气校正简化模型(SMAC),采用半经验公式代替复杂的物理过程,兼具经验方法和物理模型的优点。

目前,国内外已有多位学者对卫星影像的大气校正方法进行了研究。Adlergolden等^[5-6]基于MODTRAN4模型分别对AVIRIS短波光谱影像及HJ-1卫星CCD数据进行大气校正,得到了较好的效果。Tachiiri等^[7-8]利用6S模型分别对NOAA/AVHRR数据和GF-1数据进行了大气校正研究,并计算了影像的归一化植被指数(NDVI)值。彭妮娜等^[9]基于同步MODIS影像和6S辐射传输模型对QuickBird-2卫星海岸带影像进行大气校正,同时考虑了混浊水域上空气溶胶光学特性和高分辨率影像临近像元效应,得到了理想的结果。Fuyi等^[10]分别利用ACTOR2、ACTOR3及FLAASH模型对Landsat 5 TM影像进行大气校正研究,结合实测值分析,得出ACTOR3精度高于另外两种模型的结论。Proud等^[11]对比了6S模型与SMAC两种大气校正模型,发现SMAC具有较高的精度。宋晓宇等^[12]基于AVIRIS数据评价了大气纠正校正模块FLAASH,认为在大气条件较好的情况下,FLAASH模型反演天底点垂直观测获取的影像地表反射率更准确。周莉等^[13]采用6S模型与近红外波段离水反射率模型相结合的方法,通过神经网络模拟大气辐射传输过程,实现了可见光波段的大气校正,结果表明该方法在太湖有较好的适用性。

本文以农作物生长过程监测为应用目标,以Sentinel-2卫星影像的大气上层表观反射率(TOA)产品为数据源,研究基于SMAC、6S模型和Sen2cor方法的Sentinel-2卫星影像大气校正,主要内容包括:Sentinel-2卫星影像概况及应用潜力分析;SMAC、6S模型和Sen2cor方法原理介绍;基于SMAC、6S模型和Sen2cor方法的Sentinel-2卫星影像大气校正结果分析及精度评价。

2 数据源

2.1 Sentinel-2卫星影像概况及应用潜力分析

Sentinel-2环境监测卫星是欧盟委员会(EC)和欧洲航天局(ESA)共同倡议的全球环境与安全监测系统“哥白尼计划”中的第二颗卫星。Sentinel-2A和Sentinel-2B分别于2015年6月和2017年3月成功发射。Sentinel-2卫星的主要任务是对全球陆地表面进行高分辨率多光谱成像。该卫星的轨道高度为786 km,观测幅宽为290 km,使用寿命为7.25 a,单星重访周期为10天,A/B双星重访周期为5天。卫星上搭载了一个采用推扫式成像模式的多光谱成像仪(MSI),且包含一个光瞳直径为150 mm的望远镜,该望远镜的结构和反射镜利用碳化硅材料制成,可减少光学系统的热变形^[14]。该卫星影像共有13个波段,空间分辨率分别为10,20,60 m,各波段空间分辨率的详细情况如表1所示。

目前能够免费获取并用于植被定量监测的卫星影像多数可覆盖4个光谱波段,即蓝、绿、红、近红外波段,但许多种植物有着自身特殊的光谱,是这4个光谱波段无法表征的^[15-16]。例如,葡萄冠层的光谱颜色很特殊,水稻植株在稻田中时常被水的颜色波段所掩盖。因此,仅有4个光谱波段无法充分反映植被冠层情况。Sentinel-2卫星影像除了具有蓝、绿、红、近红外这4个光谱波段外,还具有3个植被红边波段(表1)。红边波段是指绿色植物在红波段与近红外波段之间反射率增加最快的波段,也是一阶导数光谱在该区间内的拐点。红边波段往往与植被覆盖度和叶面积指数密切相关,植被覆盖度越高,叶面积指数越大,红边斜率也就越大,相应的植被生长状态越好,红边位置会出现红移;反之,则红边位置会出现蓝移^[17]。红边波段还对植被冠层的叶绿

素和氮含量非常敏感,这对于植被定量遥感监测具有重要意义。此外,Sentinel-2卫星影像与Landsat 8 OLI、SPOT5有相似波段,在有自己独特波段和分辨率的基础上,结合其他卫星影像,进一步提高了该卫星影像的可用性。

Sentinel-2卫星提供的数据产品为Level-1C,它是经过辐射校正和几何校正的大气上层表观反射率,每幅Level-1C产品由一景100 km²的正射影像(UTM/WGS84)构成,影像的地图坐标经过数字高程模型(DEM)进行纠正,包括陆地、水、云掩模等数据。目前可以从官方网站(https:∥scihub.copernicus.eu/dhus/#/home)免费下载Level-1C数据产品。

2.2 研究区概况和数据获取

Sentinel-2卫星影像获取时间为2017年3月29日11:16,影像地理范围为115°30'E—117°30'E和39°30'N—40°00'N,具体地理位置如图1(a)所示。该影像位于北京市,影像获取时天气晴朗,该地理范围内的地形主要为平原和山地,土地利用类型有农田、林地、居民点、道路和水域。图1(b)~(d)分别为野外实验地点不同分辨率的局部放大图,图1(b)为真彩色合成,图1(c)、(d)为假彩色合成。由图可见,10 m分辨率的影像地物更清晰,内容也更加丰富;20 m分辨率的影像主要是红边波段,对于植被的长势参数非常敏感,在结合10 m分辨率影像的基础上有广阔的应用前景;60 m分辨率影像较前两者而言分辨率较低,主要用于提供气溶胶、水体和云等方面的信息。

野外实测数据采集时间为2017年3月31日11:00—12:00,天气情况与2017年3月29日基本一致。使用的探测仪器为美国Spectra Vista生产的SVC HR1024i便携式地物光谱仪,其探测波长范围为340~2510 nm,可实时获取地物辐照度及反射率曲线。SVC HR1024i便携式地物光谱仪镜头视场角为4°,当在离地面1.5 m的位置进行测量时,可视地物的面积约为0.008 m²,认为光谱仪的观测对象为单一地物。实验地点在中国农业大学上庄实验站,实验对象包括4种地物:冬小麦、裸土、水泥路和油菜地,每两个点之间的距离大于15 m,共计39个点,其分布情况如图1所示。2017年3月31日,北方植物刚开始返青,冬小麦与油菜地的植被覆盖度不高,此时卫星影像的像元为植被与土壤的混合像元。为了解决混合像元问题,在测量地物光谱的同时,利用相机拍摄了该点的照片,并裁剪出0.5 m×0.5 m的覆盖范围。利用判断条件(绿波段的DN值大于红波段和蓝波段)将照片分为植被和土壤,如图2所示,通过计算植被覆盖度,得到混合像元中植被与土壤所占的比例。在冬小麦和油菜地中,植株的种植密度及长势基本一致,可以认为该影像像元大小内的植被长势均匀,即通过照片计算的覆盖度与像元范围内的覆盖度一致。经过计算,冬小麦的植被覆盖度为0.5~0.65,油菜地的植被覆盖度为0.4左右。

图 1. 研究区位置及影像概况。(a)研究区位置;(b)分辨率为10 m,RGB通道为波段2,3,4;(c)分辨率为20 m,RGB通道为波段5,7,12;(d)分辨率为60 m,RGB通道为波段1,9,10

Fig. 1. Location and image situation of study area. (a) Location of study area; (b) resolution is 10 m, and RGB channels are band 2, 3 and 4; (c) resolution is 20 m, and RGB channels are band 5, 7 and 12; (d) resolution is 60 m, and RGB channels are band 1, 9 and 10

下载图片 查看所有图片

图 2. 图像分割。(a)实地拍摄照片;(b)图像分割结果(植被为白色,土壤为黑色)

Fig. 2. Image segmentation. (a) Photo taken in the fieldwork; (b) image segmentation result (vegetation is white, and soil is black)

下载图片 查看所有图片

3 大气校正方法

3.1 SMAC原理

大气校正简化模型(SMAC)是Rahman和Dedieu^[18]提出的,其以5S(Simulation of Satellite Signal in the Solar Spectrum)模型为基础,采用各个波段的系数方程及半经验模型来描述太阳辐射在穿越大气过程中大气介质的吸收和散射作用。SMAC的主要优点包括:1)不依赖于查找表,计算速度快;2)模型易迁移到其他辐射传输模型中;3)有正向和逆向两个过程,即大气上层表观反射率(TOA)与大气下层表观反射率(BOA)可以相互计算。当前SMAC支持的卫星包括Sentinel-2、Landsat系列MODIS、SPOT等20余种常用卫星影像。该模型同时考虑到大气介质的吸收和散射以及表观反射率的空间关系。在给定太阳天顶角、卫星天顶角、太阳-卫星相对方位角、传感器光谱响应函数和大气特性的情况下,SMAC可以实现TOA与BOA的相互计算。上述参数中,太阳和卫星的几何关系如图3所示,其中θ_s为太阳天顶角,θ_v为卫星天顶角,n为地面法线,Δϕ为太阳-卫星的相对方位角。

假设大气下层表观反射率为ρ_BOA的目标物体为均一的朗伯表面,在标准大气压下,则卫星接收到目标物的大气上层表观反射率可表示为

图 3. 太阳与卫星的几何关系

Fig. 3. Geometric relationship between sun and satellite

下载图片 查看所有图片

ρTOA(θs,θv,Δϕ)=tg(θs,θv){ρa(θs,θv,Δϕ)+[exp(-τ/μs)+td(θs)]×ρBOAexp(-τ/μv)+ρetd(θv)1-ρeS},(1)

式中μ_s=cos θ_s,μ_v=cos θ_v,t_g为大气介质总透射率(向上和向下),ρ_a为大气反射率,τ为大气光学厚度,exp(-τ/μ_s)、exp(-τ/μ_v)为直接大气透射率,td(θ_s)、td(θ_s)为大气散射透射率,ρ_e为目标物的背景反射率,S为大气球形反照率,1-ρ_eS考虑了陆地表面与大气之间的多次散射。从严格意义上来说,该方程只对单波段光谱有效,为了表达清晰,省略了波长信息^[19]。

为了降低算法的复杂性并提高计算效率,(1)式可简化为

ρBOA(θs,θv,Δϕ)=ρ'stg(θs,θv)T(θs)T(θv)+Sρ's,(2)

式中T(θ_s)、T(θ_v)表示向下或向上的整层大气透射率,ρ'_s可表示为

ρ's=ρTOA(θs,θv,Δϕ)-tg(θs,θv)ρa(θs,θv,Δϕ)。(3)

在大气传输情况下,大气介质总透射率t_g可表示为

tg(θs,θv)=exp[a(mU)n],(4)

式中m=1/μ_s+1/μ_v,U为大气垂直积分吸收量,a、n为根据6S模型数据分析特定光谱带得到的系数,这些系数可以从存储的系数表中查到。大气透射率T(θ_s)(θ=θ_s‖θ_v)同样可以通过预设系数进行简化,可表示为

T(θs)=a0+a1τ550μs+a21+μs,(5)

式中a₀、a₁、a₂为根据6S模型数据分析特定光谱带得到的系数,τ₅₅₀为波长在550 nm处的气溶胶光学厚度^[10]。

SMAC输入参数包括系数及太阳与卫星的几何参数,系数可以从SAMC的系数表中获取,影像的太阳天顶角为34.39°,太阳方位角为155.01°,各波段的卫星天顶角和卫星方位角可从影像的头文件中获取,如表2所示。

当参数输入后,经过SMAC计算得到三个系数x_a、x_b和x_c,从而得到地表反射率ρ_BOA,即

y=ρTOA-xa,(6)ρBOA=yxb+yxc。(7)

3.2 6S模型原理

6S模型是5S模型的改进版本,加入了最新的散射计算方法,太阳光谱波段的散射计算精度有所提高。6S模型可以模拟平面观测,考虑到较高的目标及非朗伯表面的边界条件,将新的气体(甲烷、一氧化二氮和一氧化碳)集成到气体传输计算当中,使得瑞利散射和气溶胶散射效应的计算精度得到显著提高^[20]。

假设地球表面为均匀的朗伯表面,在大气垂直均匀变化的条件下,基于6S模型卫星传感器测得的大气上层表观反射率定义为

ρTOA=πLλd2μsEs,(8)

式中L_λ为传感器观测到的随波长变化的大气上层辐亮度,具有方向性,与方位角、观测角有关;d为日地相对距离,可以根据日期查找;E_s为大气上层的太阳辐射通量密度,可以从Sentinel-2卫星影像的头文件中获取;μ_s为太阳天顶角的余弦值。

6S模型的参数除了包括太阳与卫星的几何参数(表2)外,还包括大气模式、气溶胶模式等参数,其中大气模式为中纬度冬季(MLW),气溶胶模式包括大陆型气溶胶模式。气溶胶厚度使用能见度模式,默认为40 km,海拔高度为0.175 km,卫星高度为786 km,光谱响应函数利用各波段的上下限及高斯函数进行模拟得到。当参数输入后,经过6S模型计算得到三个将表观辐亮度值转换为地表反射率的转换系数x_a、x_b、x_c,再通过L_λ计算得到ρ_BOA,即

y=Lλxa-xb,(9)ρBOA=y1+yxc。(10)

3.3 Sen2cor方法原理

Sen2cor是欧洲航天局(ESA)官方提供的Sentinel-2卫星影像的大气校正方法^[21],其目的是将大气上层表观反射率转化为大气下层地表反射率。该算法的核心是大气辐射传输模型libRadtran^[22]。libRadtran是用于计算地球大气辐射传输的软件包,该算法运用离散纵标法计算地球大气、太阳以及热辐射传输,包含平面、二流、多流、伪球面和球面等不同计算模块^[23],适用于各种大气条件,包括太阳和部分地球光谱的辐射、辐照度和光化通量。辐射方程可以表示为

dLβds=-L+J,(11)

式中L为位置(x,y,z)的辐亮度;β为体积消光系数;J为源函数,可定义为

J=ω4π∫p(Ω,Ω')LdΩ'+(1-ω)B(T),(12)

式中ω为单次散射反照率,p(Ω,Ω')为从方向Ω'到Ω的相位函数,B(T)为普朗克函数。当波长小于400 nm时,通常忽略普朗克常数。

Sen2cor将libRadtran代码生成的大气校正查找表集成在程序中,缩短了程序运行的时间并简化了操作流程。该查找表包含了各种大气条件、太阳角度、地形、4种基本气溶胶类型(农村、城市、沙漠和海洋)以及6种不同的臭氧含量等,存储的反射率的光谱分辨率为0.6 nm,随后利用滤波函数对光谱进行重采样,得到特定通道的反射率,查找表结构如表3所示。

三种模型的基本原理都是辐射传输理论,但SMAC基于经验的简化5S模型,6S模型则严格遵循辐射传输模型。相对于5S模型,6S模型考虑了地面目标的海拔高度、非朗伯平面的情况和新的吸收气体种类,虽然运算速度低于SMAC,但精度较高。Sen2cor基于libRadtran代码生成的查找表进行大气校正,并利用Python编译的可移植代码,提供用户命令行和集成软件两种模式,与SMAC和6S模型相比,其优点是普适性好,不需要人为输入复杂的参数,操作简单,缺点是需占用计算机的大量内存,计算速度慢。

4 结果分析

为了初步分析三种大气校正方法的结果,选取10 m分辨率(B2,B3,B4,B8)及20 m分辨率(B5,B6,B7,B8a)8个波段,并选择冬小麦、裸土、水泥路和油菜地4种地物的像元实测值、大气上层表观反射率及三种模型的校正结果进行分析。其中冬小麦和油菜地为混合像元,首先通过第2.2节的方法计算植被覆盖度,得到混合像元中植被与土壤的比例,并结合实测地物光谱计算出混合光谱,结果如图4所示。

由图4可见,Sentinel-2卫星影像经过三种大气模型校正后,其光谱曲线与实测地物的光谱曲线的变化趋势一致,有较高的拟合度,蓝波段Band 2经大气校正后相比于大气上层表观反射率TOA明显下降,红边、近红外波段明显上升。在冬小麦和油菜地的光谱曲线中, Band 2的三种大气校正结果高于实测地物的反射率,误差较大。Sen2cor方法与SMAC的误差约为0.04,6S模型的误差约为0.06;

图 4. 各类地物大气校正前后的反射光谱曲线。(a)冬小麦;(b)油菜地;(c)裸土;(d)水泥路面

Fig. 4. Reflectance spectral curves of different ground objects before and after atmospheric correction. (a) Winter wheat; (b) rape field; (c) bare soil; (d) cement pavement

下载图片 查看所有图片

Band 3~Band 5的大气校正结果与实测值拟合效果较好,误差均小于0.02;在大于705 nm的Band 6~Band 8a波段中,实测地物的反射率高于三种大气校正结果,误差范围为0.01~0.04。在土壤光谱曲线中,大气校正结果相比于大气上层表观Band 2反射率明显下降,Band 4~Band 8a波段的反射率明显增加,Band 2的三种大气校正结果高于实测地物的反射率,其余波段的实测地物反射率均高于三种大气校正结果,误差在0.02左右。在水泥路的光谱曲线中,大气校正结果相比于大气上层表观反射率除了Band 2明显下降,其余波段反射率明显增加,而且实测地物的反射率均高于大气校正结果。从图4(d)中可以看出,绿波段Band 3、红波段Band 4和红边波段Band 5的误差较大,三种大气校正结果的误差范围为0.04~0.07,其余波段的误差较小,误差范围为0.02~0.04。图4中的结果表明,4种地物大气校正的精度从高到低依次为土壤、植被(冬小麦和油菜地)、水泥路,其误差可能源于混合像元、气溶胶类型和测量误差。

SMAC、6S模型和Sen2cor方法都是基于大气辐射传输理论的机理模型,只是计算方法不同。为了验证三种模型的一致性,对其结果分别进行回归分析,图5和表4分别为三种大气校正结果之间的相关性分析和回归分析。由图5可见,用三种方法校正同一影像时,结果具有很强的相关性。由表4可见,6S模型与Sen2cor方法的决定系数R²最高(0.997),均方根误差E_rmse最小(0.0031);SMAC与Sen2cor方法的决定系数R²为0.9906,均方根误差E_rmse为0.006;6S模型与SMAC的决定系数R²为0.9874,均方根误差E_rmse为0.0063。三种结果的检验P值均小于0.0001,可信度在99.99%以上。结果表明三种模型具有较好的一致性。

图 5. 基于Sen2cor,6S模型和SMAC的大气校正结果相关性分析。(a) SMAC, 6S与Sen2cor的反射率结果; (b) SMAC与Sen2cor的反射率结果; (c) SMAC与6S的反射率结果; (d) 6S与Sen2cor的反射率结果

Fig. 5. Correlation analysis of atmospheric correction results of Sen2cor, 6S model and SMAC. (a) Result of reflectance by SMAC, 6S and Sen2cor; (b) result of reflectance by SMAC and Sen2cor; (c) result of reflectance by SMAC and 6S; (d) result of reflectance by 6S and Sen2cor

下载图片 查看所有图片

为了定量分析SMAC、6S模型和Sen2cor方法大气校正的精度,将野外39个实测地物的光谱数据与大气校正结果提取出Band 2~Band 8a范围的反射率进行回归分析,结果如图6和表5所示。由图6可见,三种模型均与实测值有较强的相关性,且数据整体较集中。当反射率为0.2时,数据相对离散,可能的原因是此时为红边波段,即光谱曲线中变化率最大的波段,误差可能较大。从拟合的趋势线可以看出,Sen2cor方法的大气校正效果优于6S模型和SMAC。从表5可以得出,Sen2cor方法与实测数据的决定系数R²最高(0.8196),均方根误差E_rmse最小(0.0388);SMAC与实测数据的决定系数R²最小(0.8137),均方根误差E_rmse为0.0395。

三种结果的检验P值均小于0.0001,可信度在99.99%以上。

图 6. 反射率大气校正结果与实测值的相关性分析

Fig. 6. Correlation analysis between atmospheric correction results and measured values of reflectance

下载图片 查看所有图片

植被定量遥感是卫星影像最普遍的应用之一。植被指数是常用的参数,其通常是由两个或两个以上的波段计算得到。大气校正的目的是让卫星影像的光谱趋近于地表真实光谱,所以植被指数也用来检测大气校正效果^[24]。本研究选用归一化植被指数(NDVI)作为指标,对大气校正后影像上的NDVI值及实测值进行分析,NDVI值可表示为

fndvi=ρNIR-ρredρNIR+ρred,(13)

式中ρ_NIR、ρ_red分别表示近红外波段和红光的反射率,结果如图7所示。其中SMAC的NDVI值最大,其次为Sen2cor方法和6S模型。表6为Sentinel-2影像NDVI与实测NDVI值的回归分析,由表可知,SMAC计算的NDVI值与实测值的相关性最高,决定系数

R²为0.6389,均方根误差E_rmse为0.093。三种结果的检验P值均小于0.0001,可信度在99.99%以上。

图 7. Sentinel-2影像NDVI与实测NDVI值的相关性分析

Fig. 7. Correlation analysis between NDVI from Sentinel-2 imagery and measured NDVI value

下载图片 查看所有图片

5 结论

在对Sentinel-2卫星影像和大气校正理论研究的基础上,基于SMAC、6S模型和Sen2cor方法对Sentinel-2卫星影像进行了大气校正,并结合实测地物光谱对其精度进行了分析。结果表明,影像光谱曲线与实测地物的光谱曲线的变化趋势一致,二者具有较高的拟合度。从三者之间的相关性分析得出,三种模型具有较好的相关性和一致性,其中6S模型与Sen2cor方法的决定系数R²最高,SMAC与Sen2cor方法次之,6S模型与SMAC最小。从Band 2~Band 8a波段的大气校正结果来看,Sen2cor方法的精度最高,其次为6S模型和SMAC。从NDVI值的分析可以看出,SMAC计算的NDVI值与实测值的相关性最高,其次为6S模型和Sen2cor方法。三种方法的大气校正效果显著,且具有较高的精度,Sentinel-2卫星影像经过校正后影像质量得到明显提高,增加了可用性。

在后续的研究中,将着重分析大气校正参数(如气溶胶类型)对其结果的影响,还将选择不同季节、不同植被类型和植被覆盖度的数据进行验证。