大气校正是进行地表参数定量分析和遥感反演前必不可少的处理步骤,其目的是为了消除大气分子、气溶胶的散射以及大气中H₂O,O₂,CO₂,CH₄和O₃等物质的吸收对地表反射率、辐亮度等相关物理量的影响。在遥感数据的应用过程中为尽量消除大气影响,常将遥感数据定标后的天顶反射率转换为反映地物真实信息的地表反射率^[1,2,3]。目前较为常见的大气校正方法有: 基于辐射传输模型的MORTRAN模型、LOWTRAN模型、ATCOR模型、6S模型和FLAASH模型,基于简化辐射传输模型的黑暗像元法,基于统计学模型的反射率反演法,基于统计的不变目标法和直方图匹配法等^[4]。其中,6S(second simulation of the satellite signal in the solar spectrum)模型是在5S模型(simulation of the satellite signal in the solar spectrum)基础上发展起来的一种对大气进行精确校正的模型^[5],采用逐次散射法 (successive orders of scattering,SOS)和最新近似法(state of the art)计算散射和吸收,综合考虑了高程、非朗伯体表面以及新的吸收物质(CH₄,NO₂,CO)等条件,提高了瑞利散射与气溶胶散射的计算精度^[6]。FLAASH(fast line-of-sight atmospheric analysis of spectral hypercubes)模型则基于MODTRAN4辐射传输模型,算法精度高,任何有关遥感图像的标准MODTRAN大气模型和气溶胶类型都可以直接使用,可有效去除水蒸气和气溶胶的散射效应。

目前国内对SPOT6图像的大气校正研究较少,且对其在高海拔地区适用的大气校正方法也没有较为深入的研究。因此,本文主要利用6S及FLAASH模型,辅以相关大气数据对覆盖研究区的SPOT6图像进行大气校正; 其中对6S模型将依据研究区的气溶胶光学厚度(aerosol optical depth,AOD)和海拔的平均值及梯度值划分为AVG6S模式和GRD6S模式2种模式对图像进行处理; 最后使用Landsat8 SR地表反射率产品作为验证数据,通过相关性分析、典型地物光谱曲线和归一化植被指数(normalized difference vegetation index,NDVI)对比等指标进行精度验证,以期得到适用于SPOT6的高海拔地区最佳大气校正方法。

湟水流域地处青藏高原与黄土高原过渡地带,地理位置介于E100°42'~103°04',N36°02'~37°28'之间(图1),流域面积约为16 120 km²。区内地形由河谷冲积平原、黄土丘陵和中高山地组成,西北高东南低,海拔在1 662~4 860 m之间。地貌类型多样,高程变化剧烈,AOD分布不均匀。根据青海省气象局大气监测站点数据,2016年8月份主要影响空气质量和能见度的PM₁₀,PM_2.5和SO₂含量分别为80.06 μg/m³,45.03 μg/m³和16.64 μg/m³,无严重大气污染; O₃,NO₂和CO的含量分别为123.09 μg/m³,32.01 μg/m³和0.75 mg/m³。

覆盖整个湟水流域约需5景SPOT6图像,因篇幅限制,只选择其中1景作为本文论述对象和研究目标。此区域位于西宁市东部,土地覆盖类型包括城市、河流、农田、草地丘陵和高山林地等(图1(b)),具有一定代表性和典型性。

SPOT6卫星于2012年9月9日发射,拥有空间分辨率为6 m的多光谱波段(蓝光、绿光、红光、近红外)和1.5 m的全色波段,标准图像覆盖范围为60 km×60 km^[7]。本文使用的SPOT6卫星图像为镶嵌成果,非标准幅宽,覆盖范围约为5 528.9 km²。与其他高空间分辨率卫星相比,SPOT6具有成像图幅宽、重访周期短等特点。

对于大气校正后的地表反射率产品最好的验证数据为光谱仪野外观测获取的实测光谱,再经计算转换模拟得到的波段等效遥感反射率^[8]。但对于历史图像数据,上述方法是不可行的。同日期的MODIS地表反射率产品MOD09通常被用来作为地表反射率产品验证数据,但对于SPOT6图像,500 m的空间分辨率显得过低,且对于地形复杂的本文研究区,可能包括大量的混合像元和复杂下垫面。因此,本文选取Landsat8 SR地表反射率产品作为验证数据^[9],其与SPOT6对应波段光谱宽度相近,且此景SPOT6图像数据的获取时间为2016年8月8日,Landsat8 SR的获取时间为2016年8月7日,天气状况基本一致。

在进行大气校正之前,需要先对SPOT6原始图像进行辐射定标预处理。辐射定标就是将图像的DN值化为辐射亮度值或表面温度等物理量的处理过程^[10],计算公式为

式中: L_i为第i波段的辐射亮度; DN_i为第i波段的灰度值; Gain_i和Bias_i为第i波段的定标系数。

大气层顶表观反射率(top of atmosphere apparent reflectance,TOA)是指大气层顶的反射率,它等于地表反射率加上大气反射率,故可理解为大气校正前的反射率。TOA的计算公式为^[11,12]

式中: L_λ为波段λ处(或第i波段)的辐射亮度; D为日地距离; E_SUNλ为波段λ处(或第i波段)的大气顶层太阳辐照度; θ为太阳天顶角。

6S模型是一个基于物理原理的模型,没有在一个特定的卫星场景、测试场地上进行过优化,只需要正确地输入参数,使其有利于对遥感数据进行比较或全局分析^[13]。所以,6S模型输入参数的合理与否直接影响着大气校正结果质量的好坏。输入参数主要包括几何参数、气溶胶模式、大气模式、光谱响应函数、目标海拔和传感器高度等^[14]。本文主要考虑气溶胶模式和目标海拔这2个因素。因为,在整个大气传输过程中,CO₂和O₂在全球范围内的变化不大,且在短波宽带VNIR传感器中作用较小,对其吸收作用可以忽略。在大气窗口区域发生辐射传递的主要大气参数变量是水蒸气含量以及气溶胶类型和AOD^[15]; 而对气溶胶散射的校正是大气校正的重点和难点,因为AOD很难获取,且在1景图像中不同地点的AOD是不同的,本景SPOT6图像中,城市上空的AOD比农林地区域高出约2倍。获取研究区AOD最好的办法是利用太阳分光光度计现场测量,对于覆盖面积广、下垫面复杂的图像,需要测量很多样点; 而对于历史图像数据,此种方法是不可行的^[16]。所以,本文使用MOD/MYD04_3K产品作为获取研究区AOD的数据源。此外本研究区海拔变化剧烈,最大高程差可达2 481 m,而O₃的吸收作用和海拔却有较高的依赖关系,在分子反射贡献已知的情况下,其反射强度取决于压力水平和海拔^[17],也就是说地表反射率在复杂地形区会受到地形条件影响^[18],这也给本次研究增加了困难。

综合考虑上述2个因素,本文采用了2种6S模型的参数输入模式。一种是整景采用平均海拔和平均AOD(下文表示为AVG6S模式); 另一种则根据AOD差和海拔差设置梯度气溶胶和梯度海拔建立查找表(下文表示为GRD6S模式)进行大气校正。表1列出2种6S模式的输入参数。

将表1中的参数作为6S模型的输入参数,对SPOT6多光谱数据的4个波段进行大气校正计算,得到大气校正参数x_a,x_b和 xc[19],对于GRD6S模式,单个波段有35组大气校正参数,每个像元再通过插值得到对应的校正参数,然后计算校正后的地表反射率,即

式中: ρ为校正后的反射率; x_a,x_b和x_c为大气校正参数; L为经过辐射定标的大气顶层辐亮度。

根据上述数据和模型参数,利用IDL语言实现对SPOT6多光谱数据的大气校正^[20],具体流程见图2。

选取城区和农村地区2块区域进行直观定性分析。由于SPOT6图像为交付产品,整体质量较好,没有明显的大气污染现象。所以,在AOD大的城区用GRD6S模式校正出来的结果更加清晰,对比度增加,居民点和城区植被接近真实颜色,效果明显(图3)。而在农村地区,3种校正结果并没有很显著的视觉对比(图4),主要原因可能是农村地区雾霾等大气污染并不严重,AOD本来就小,加之图像质量良好。

3.2.1 校正结果与Landsat8 SR的线性回归分析

为了验证校正后SPOT6图像的质量,将Landsat8 SR与TOA,GRD6S模式,AVG6S模式和FLAASH模型的处理后图像进行对比,随机选取200个像素点绘制散点图,做线性回归分析,以此对比两者间的差异,并计算出决策系数R²和均方根误差(root mean square error,RMSE)(图5)。

从图5可以看出,TOA,GRD6S模式,AVG6S模式和FLAASH模型与Landsat8 SR的R²分别为0.955,0.977,0.976和0.974,3种大气校正结果均与Landsat8 SR有较强的相关性; RMSE分别为0.018,0.004,0.002和0.006,均小于0.01,说明3种大气校正方法均能较好地校正SPOT6原始图像,校正后的图像质量更好。

3.2.2 不同校正方法处理前后分析

为了更好地分析3种校正方法的优势和差异,本文根据地形要素和大气状况分别选取了3块试验区: ①海拔大于3 000 m,AOD小于0.2的高山区; ②海拔在2 600~2 900 m,AOD在0.3左右的农村地区; ③海拔小于2 300 m,AOD大于0.45的城区。三者的面积分别为498 km²,196 km²和515 km²。表2统计了3块试验区的TOA和3种大气校正方法结果的标准差(Std)和其与Landsat8 SR的相对误差(relative error,RE),如表2所示。标准差反映了图像离散程度,标准差越大,灰度级分布越分散,图像的目视效果越好; RE反映了校正结果的可信度,其值越小,则校正质量越好,RE定义为

式中: V₁为TOA以及3种校正结果的反射率; V₂为Landsat8 SR的反射率。

由表2可知,整体而言,整景图像3种校正结果和Landsat8 SR的4个波段的标准差都大于TOA的标准差,说明校正后图像的区分度比校正前的好,图像信息得到增强; 3种校正结果的4个波段的RE均远小于TOA对应波段的RE,且随波长的增加总体呈现递减趋势,主要是因为瑞利散射原理,波长越短受大气的影响越大,随着波长的增加受大气的影响逐渐变小。

1)对整景图像而言,蓝、绿波段,AVG6S模式的RE最小,分别为2.79%和7.97%,表现最好; 3种校正结果在红波段的RE分别为: 0.92%,1.62%和12.29%,GRD6S模式和AVG6S模式表现优于FLAASH; 在近红外波段并无太大差异,RE分别为3.01%,1.91%和1.15%。因此从整景图像上看,GRD6S模式和AVG6S模式的校正效果整体上要好于FLAASH模型,表现最好的是AVG6S模式。

2)对各试验区而言,城区低海拔高AOD,所以校正前后对大气状况敏感的蓝波段对比效果最为明显,GRD6S模式的RE远低于其他两者,为4.15%,对蓝波段的处理最好; 其他3个波段,6S模型的2种模式也整体好于FLAASH模型。农村地区AVG6S模式表现最好,其次为GRD6S模式和FLAASH模型。高山区FLAASH模型在蓝、绿波段表现更好,RE分别为23.93%和8.66%,其次是GRD6S模式,而AVG6S模式的RE最大,分别为37.38%和18.84%,其原因可能是高山区各像元的海拔和AOD与AVG6S模式所设置的海拔和AOD参数相差较大,从而导致在本区域内与Landsat8 SR相差较大; 而在红、近红外波段GRD6S模式表现更好,其RE分别为1.25%和0.73%。

综上所述,6S模型在研究区的校正效果要好于FLAASH模型。而对于GRD6S模式和AVG6S模式而言,分别有各自的优势。从整景看,AVG6S模式的校正效果是最好的,主要原因是与所设平均AOD和平均海拔相似的区域较广,约占总面积的72.6%; 而在城区和高山区,综合考虑,GRD6S模式的校正效果更好些。FLAASH模型逊于6S模型的原因可能是,前者更多地使用图像本身的信息来进行大气校正,而对于本研究区复杂的地理环境,使用图像当天大气状况参数建立的6S模型,能进行更准确的大气校正。

3.3.1 植被

健康植物的光谱曲线有明显的特点,在可见光的0.55 μm附近有一个反射率为10%~20%的小反射峰,在0.45 μm和0.65 μm附近有2个明显的吸收谷,在0.7~0.8 μm有一个陡坡,反射率急剧增高,并在近红外波段形成一个反射峰^[21,22]。为研究3种校正方法的特征及差异性,依旧针对高山区、农村地区和城区选取植被进行对比,结果如图6所示。光谱曲线结果由随机选取60个采样点的平均值所得。

从图6可以看出,校正前TOA的蓝波段反射率略高于绿波段,并没有表现出上文所描述的“峰”和“谷”,也不符合植物的典型地表反射率规律; 而校正后图像的光谱曲线则表现出典型的植物光谱特征。高山区FLAASH模型在蓝绿波段好于6S模型,而在红波段和近红外波段不如6S模型; 在农村地区6S模型明显好于FLAASH模型,尤其是AVG6S模式的光谱曲线与Landsat8 SR吻合度最高,GRD6S模式也基本吻合; 而在城区GRD6S模式和AVG6S模式也优于FLAASH模型,GRD6S模式好于AVG6S模式,虽然在绿波段不如AVG6S模式,但在其他3个波段更接近验证波段,尤其是受气溶胶影响较大的蓝波段。所以从植被光谱对比可以看出,无论是AVG6S模式还是GRD6S模式,校正效果都要好于FLAASH模型; 对于前两者,AVG6S模式在农村地区表现最好,而GRD6S模式则在城区表现更为出色,其主要原因是城区AOD大于全景平均AOD,海拔也低于平均海拔,GRD6S模式所采用的梯度参数更接近于城区实际情况,GRD6S模式的优势得以体现。

3.3.2 土壤

典型土壤的反射率光谱特征一般随波长增加而变大; 颜色较浅、土壤颗粒细微的土壤反射率较高,且土壤表面覆盖度小于15%时,其光谱反射体征与裸土相近^[21,22]。本文的研究区在青藏高原与黄土高原的过渡带,东部与黄土高原接壤,土壤特征与上述表述相近。图7由100个土壤采样的地表反射率平均值所得。可以发现,3种校正方法所得土壤光谱反射率曲线总体趋势与验证光谱曲线一致,均能很好地反映真实土壤光谱情况。由于黄土区多位于低山丘陵区,海拔和AOD均接近平均值,所以AVG6S模式的校正效果较好,其次为GRD6S模式、FLAASH模型。

3.3.3 水体

不同水体因其各异的水体特性会造成反射光谱特征的差异。在可见光范围内,水体的反射率总体上比较低,并随着波长的增大逐渐降低,波长大于0.75 μm时,水体几乎成为全吸收体; 而泥沙较多的浑浊水体,反射波谱曲线整体偏高,且随着悬浮泥沙度的增加,波谱反射峰向长波方向移动(即“红移”)^[21,22]。本文分别从水库和湟水选取采样点,水库水质较清澈,有机叶绿素含量较高; 湟水为黄河的重要支流,水体浑浊,无机泥沙含量较高。图8为大气校正前后60个采样点在各波段的平均值。校正后的水体符合上述描述,在水质较清澈的水库,除因有机叶绿素的缘故,在绿波段有一个小反射峰外,反射率随光谱波长增加而降低; 在泥沙含量较高的湟水,反射率总体较水库高且随波长的增加呈上升趋势。

从图8可以看出,FLAASH模型校正的结果更接近验证数据,而6S模型所得反射率较验证结果均偏低。说明FLAASH模型对于水体的校正效果好于6S模型。

3.3.4 城区建筑

对于城区环境,选取建筑物和公路2种地物对3种大气校正结果进行对比(图9)。建筑物屋顶和公路多为水泥材质,且灰尘较多,因而反射率较高,光谱曲线呈现整体上升趋势,在可见光波段内上升较慢,进入近红外波段反射率上升较快。从图9可以看出,校正后的光谱曲线与验证光谱趋势大体一致。三者相比较,表现较好的为GRD6S模式,其次为AVG6S模式,最后是FLAASH模型,但都明显好于TOA的光谱曲线。GRD6S模式有较好的校正效果的原因可能是所采样点在河谷城区和乡镇居多,而在此处的AOD偏大。

植被定量遥感是卫星图像最普遍的应用之一,植被指数是常用的参数。大气校正的目的是让图像的光谱接近地表真实光谱,所以植被指数也常被用来检测大气校正效果。本文选用NDVI作为检验指标。

表3为TOA,GRD6S模式,AVG6S模式,FLAASH模型与Landsat8 SR在NDVI的结果分析表,分别统计了其各自的标准差(Std)以及与Landsat8 SR的RE。

从整景图像来看,3种校正结果的RE均远小于TOA的16.86%,其中AVG6S模式最小,其次为GRD6S模式,最后是FLAASH模型。三者均得到了很好的校正效果,且2种6S校正模式要优于FLAASH模型,其中AVG6S模式校正效果最好。另外,同样分析了3块试验区各大气校正方法的校正结果。由表3可以看出,在高山区和城区GRD6S模式拥有较小的RE,校正效果最好,AVG6S模式其次。农村地区AVG6S模式相对误差较小,校正效果最好; GRD6S模式次之。由此可得,6S模型2种模式的大气校正方法对于本研究区域要优于FLAASH模型。而对于GRD6S模式和AVG6S模式2种校正模式而言,校正差异并不大,但各自有校正优势区域,结论与3.2.2一致。

在整景NDVI中随机采样60个点,绘制出TOA,GRD6S模式,AVG6S模式,FLAASH模型和Landsat8 SR的折线图(图10)。

从图10可以看出,TOA的NDVI值最小即折线在最下方,FLAASH模型其次。受采样点的分布、不同的地形和大气状况条件影响,GRD6S模式和AVG6S模式在Landsat8 SR一定范围内波动。因为大气的瑞利散射和气溶胶的散射作用在可见光波段比O₃的吸收作用强,且波长越短受影响越大,而近红外波段受强烈的水汽(吸收)效应要大于可见光波段^[23]。大气校正削弱了大气瑞利散射和气溶胶散射而导致的可见光波段信号的增强,并补偿了因为大气衰减而引起的近红外波段减小的信号^[24]。两波段的对比度较校正前增高,图像的NDVI也随之增高。3种大气校正方法的结果也都符合上述分析。较之FLAASH模型,2种6S模式得到更宽的NDVI值域,尤其NDVI在0.8以上的区域多为本文研究区中特有的青海云杉林区和高寒草甸区,校正后增加了区分度,有利于后期的植被分类研究,是更适合本区域的大气校正方法。

本文以青藏高原和黄土高原过渡区的湟水流域为研究区,采用6S模型的2种模式和FLAASH模型对SPOT6图像进行大气校正,以获取日期和天气状况相近的Landsat8 SR为验证数据,进行对比分析研究。得到如下结论:

1)6S模型和FLAASH模型均能消除大气因素的影响,得到较为准确的地表反射率结果。相对而言,在水体的校正上FLAASH模型效果好于6S模型,在其他典型地物的反射率光谱曲线对比和NDVI检验中,6S模型整体上都优于FLAASH模型,能较好反映地物的真实地表反射率,说明6S模型是最适宜于本研究区的大气校正方法。

2)对于6S模型的2种模式,从整景影像分析,AVG6S模式的校正效果最好。但两者在特定区域各具优势,在农村地区,AVG6S模式的4个波段RE最小; 而在城区和高山区,GRD6S模式则略好于AVG6S模式; 在NDVI对比结果分析中,其RE分别为0.23%和1.49%。所以当图像整体质量较好、AOD波动不大时,AVG6S模式校正效果最佳; 而当图像大气状况较差,AOD波动较大时,则宜采用GRD6S模式,以得到更好的校正效果。

3)通过本文的研究,对于高海拔复杂地形区的图像处理,给出以下建议: FLAASH模型不依赖遥感成像时的大气参数数据,且操作简单,有成熟商业软件可以使用,在缺少大气参数或数据不理想的时候推荐使用; 当有大气参数数据且图像本身质量较好时,推荐使用AVG6S模式校正方法; 若图像大气状况较差,则建议使用GRD6S模式校正方法; 另外,6S模型校正方法在批量处理方面优势明显。

4)本文为简化流程,仅考虑了均一地表; 且GRD6S模式只采用了简单的线性插值。这2方面的问题需要在今后研究中进一步改善解决。