地表反射率是地表物体反射能量与到达地表物体的入射能量的比值，它是反演很多地表参数的基础变量，对于研究地气间辐射能量平衡、地物识别和分类、水分与气候模型等具有十分重要的意义^［1-2］。地表反射率的计算对于植被覆盖变化监测和地表物理参数反演显得尤其重要，包括叶面积指数、叶绿素含量和其他生物量反演等。遥感影像大气校正的目的是去除大气分子、气溶胶和水蒸汽等大气成份吸收和散射对传感器接收信号的影响，大气校正是遥感地表反射率反演中的关键环节。大气校正的研究始于20世纪70年代，经过30多年的发展，国内外学者提出了很多大气校正方法，主要包括辐射传输模型（Radiative Transfer Model）法^［3-5］、暗目标减法^［6-8］（Dark Object Subtraction）、不变目标法（Invariant Object Method） ^［9-10］、直方图匹配法（Histogram Matching Method）^［11］等。辐射传输模型法是较常用的大气校正方法，它基于电磁波在大气中的辐射传输原理，具有严密的理论基础和较高的校正精度，其算法在原理上基本相同，差异在于不同的假设条件和适用的范围，因此产生很多可选择的大气较正模型，最常用的如6S（Second Simulation of the Satellite Signal in the Solar Spectrum）模型、MODTRAN （Moderate Resolution Transmission）模型、LOWTRAN（Low Resolution Transmission）模型、空间自适应快速大气校正模型ATCOR（A Spatially-Adaptive Fast Atmospheric Correction）等。6S模型是在法国大气光学实验室Tanre等和美国马里兰大学地理系Vermote在5S（Simulation of the Satellite Signal in the Solar Spectrum）模型的基础上发展起来的^［12］。该模型是用FORTRAN语言编写的适用于太阳反射波段（0.4~2.5 μm）的大气辐射传输模型。6S模型假定晴空无云的条件下，考虑了H₂O、O₃、O₂、CO₂、CH₄、N₂O等气体的吸收、分子和气溶胶的散射以及非均一地面和双向反射的问题。6S模型采用了最近近似（state of the art）算法和逐次散射（successive orders of scattering）算法来计算散射和吸收，改进了模型的参数输入，提高了计算精度。6S模型具有较高的精度且不受研究区特点以及目标类型等的影响，应用范围很广泛。

陆地卫星生态系统干扰自适应处理系统（Landsat Ecosystem Disturbance Adaptive Processing System， LEDAPS）由美国航空航天局（NASA）Making Earth System Data Records for Use in Research Environments （MEaSUREs）项目资助开发，该系统实现了Landsat长时间系列卫星（Landsat 4、5、7）地表反射率的自动化生产。LEDAPS基于6S辐射传输模型，地表反射率产品生产所需要的大气水蒸汽（Water Vapor， WV）、气压（Surface Pressure， SP）来自于NCEP（NOAA National Centers for Environmental Prediction）再分析数据，分辨率为2.5°×2.5°，臭氧数据来自TOMS （Total Ozone Mapping-Spectrometer）数据，分辨率为1.25°（经度）×1.00°（纬度），气溶胶光学厚度（Aerosol Optical Thickness，AOT）数据则是直接根据Landsat影像采用暗像元浓密植被法（Dense Dark Vegetation， DDV）进行反演得到。目前国际主流的Landsat 4/5/7卫星地表反射率产品大多基于成熟的LEDAPS软件进行生产。Landsat 8作为一种新型卫星，其携带的OLI（Operational Land Imager）传感器在波段的数量、波段的光谱范围和影像的辐射分辨率上比Landsat 4、5、7卫星都有所改进，研发针对Landsat 8卫星的地表反射率反演算法是目前国内外的一个研究热点。USGS发布的Landsat 8地表反射率产品基于内部算法（internal algorithm）来进行大气校正。Wang等^［13］提出了一种基于MODIS数据的逐像元大气校正算法，用于生产2013年覆盖整个中国的Landsat-8 OLI 地表反射率，该算法是基于6S模型建立的，结果表明基于6S模型的逐像元大气校正是可靠的。6S辐射传输模型一般可直接进行大气校正或预先为大气校正建立查找表。直接法比较简单，但非常费时。如果将1 km空间分辨率MODIS大气产品作为输入参数来对Landsat 8 OLI影像进行逐像元大气校正，则6S模型将会被调用2.4×10⁵次，直接法不适用于大规模Landsat 8影像的快速大气校正。

为实现Landsat 8影像的快速大气校正，本文基于查找表开发6S模型逐像元快速大气校正算法。应用该方法在复杂地形山区计算地表反射率及NDVI，并与新近提出的阴影消除植被指数SEVI进行抗地形阴影影响等性能比较，对比验证6S模型逐像元快速大气校正结果在复杂地形山区植被信息反演精度。

研究区选择北京八达岭国家森林公园。八达岭国家森林公园位于万里长城八达岭和居庸关之间，总面积2 933 hm²，最高峰海拔1 238 m，分布植物约539种、动物约158种、林木覆盖率达96%，为中国首家通过FSC国际认证的生态公益林区。

研究区的Landsat 8 OLI卫星数据来自于中国遥感卫星地面站，影像的坐标系为WGS84，投影为横轴墨卡托投影，过境时间为2016年5月20日，行号为123，列号为32，空间分辨率为30 m。影像获取时的太阳方位角和高度角分别为156.793 278°和43.279 056°，选取488×541个像元子集作为研究区，该研究区的地形高程变化范围为82~1 029 m。

在6S模型的输入参数中，考虑了几何条件（包括太阳天顶角、方位角、观测天顶角和方位角）、大气成分（水汽、臭氧、气溶胶类型和浓度）、光谱波段、地面反射类型和光谱的变化等。然而，在实际使用中很难收集到实时的大气剖面数据，特别是对于历史遥感数据。为尽可能准确地获取6S模型需要的大气参数，将与Landsat 8卫星过境时间同步的Terra MODIS大气参数产品作为6S模型的输入，以提高地表反射率反演的精度。其中大气水蒸汽含量采用MODIS MOD05产品，臭氧数据来自MODIS MOD07 Total Ozone 产品，气溶胶光学厚度来自MODIS MOD04产品。总体大气校正流程图如图2所示。

Landsat 8 OLI遥感数据在反演地表反射率之前，需要对影像的DN值进行辐射定标，即根据公式（1）将影像像元亮度值转换为星上辐射亮度：

其中：Lrad为辐射亮度，Mρ为各波段的增益，从头文件获得（RADIANCE_MULT_BAND_i， i为波段号）；Aρ为各波段的偏置，从头文件获得（RADIANCE_ADD_BAND_i， i为波段号）；DN为各波段的像元值。同时，影像的大气表观反射率可以从公式（2）得到：

其中：Mρ'可从头文件获得（REFLECTANCE_MU⁃LT_BAND_i， i为波段号）；Aρ'从头文件获得（REFLECTANCE_ADD_BAND_i， i为波段号）；DN为各波段的像元值；θ为太阳天顶角。

中分辨率成像光谱仪（MODerate-resolution Imaging Spectroradiometer） -MODIS是Terra和Aqua卫星上搭载的主要传感器之一，Terra卫星每日地方时上午10：30过境，这与Landsat 8 OLI过境时间相差不到半个小时。在这项研究中，MOD04、MOD05和MOD07大气参数产品都可以从MODIS网站免费下载，分别用于提供550 nm的气溶胶光学厚度、大气水蒸汽含量和臭氧含量。由于云以及其他因素的影响，这些产品有时会存在无效的像元值，本文用Kriging插值法进行插值。

为了提高大气校正的速度，建立一个Landsat 8大气参数查找表，该查找表的设计主要涉及维度以及各维度的节点数。如果维数或节点的数目太大，查找表的建立将会非常耗时，还有可能存在数据冗余的问题。相反，当维数或节点的数量太少，计算精度会由于欠采样而受到影响。Landsat 8 OLI大气校正主要影响因素包括太阳天顶角、水汽、臭氧、气溶胶光学厚度（AOT）和目标高度。

为了得到更合理的维度和节点数，对上述因素进行了敏感性分析研究，结果如图3所示。图3中Y轴是地表反射率与初始的变化百分比。由于地面假设为朗伯面，太阳方位角对地表反射率没有影响，因此，太阳方位角被设置为固定值，其他影响因素作为查找表的维度。考虑到这些参数的分布范围，并参考Peng等^［14］的研究成果，确定5个输入参数的节点如表1所示。

根据表1，参数组合大约有386 100组，每一套组合的模拟都是基于6S V1.0B模型，为了节省空间和内存，所有组合的模拟都以二进制形式存储。利用处理后的MODIS大气参数影像、DEM数据以及Landsat 8数据的成像几何参数，基于查找表逐像元插值计算该点的大气校正系数，最后根据大气校正系数得到其地表反射率图像。其中所采用的插值算法为多维线性插值法。

根据图2所示流程，建立了基于6S辐射传输模型的Landsat 8 OLI逐像元快速大气校正算法，通过该算法，所有的Landsat 8 OLI影像都可以自动快速地转换为地表反射率影像。图4为2016年5月20日中国八达岭国家森林公园Landsat 8 OLI影像大气校正结果。由于消除了大气的影响，地表反射率影像比原始影像更加清晰。

地表反射率结果的实用性很大程度上取决于质量和准确性。通常有许多方法来进行算法检验，如典型绿色植被的光谱分析、实地测量验证与公开发布的地表反射率产品比较分析（交叉验证）等。光谱分析是一种间接的验证方法，通过研究典型对象（如绿色植被）的光谱，并利用光谱的光谱变异率（如归一化差值植被指数（NDVI））来间接证明结果的合理性。实地测量验证通常计算均方根偏差（RMSD），这是最直接的验证方式，然而，实地观测有许多局限性，如受天气和地形的限制，对于历史遥感数据，往往无法获取到相应的实测反射率数据来进行精度验证。本文利用光谱分析和交叉验证相结合来评估模型的质量，通过比较分析表观反射率、模型获得的地表反射率以及USGS产品的植被光谱以及NDVI指数来对算法进行总体评价。

植物的光谱特征可使其在遥感影像上有效地与其他地物相区分。典型植被在可见光的0.55 μm附近有一个反射率为10%~20%的小反射峰，在0.45 μm和0.65 μm附近有两个明显的吸收谷。在0.7~0.8 μm反射率急剧增高，在近红外波段0.8~1.3 μm之间形成一个高的，反射率可达40%或更大的反射峰。由于路径辐射，可见光波段的表观反射率往往高于地表反射率。在近红外和短波红外波段，由于大气的吸收，表观反射率低于真实的地表反射率。红光和近红外波段以不同的方式受到大气的影响，光谱指数NDVI由近红外波段和红光波段的比值计算得到，因此，大气校正后的NDVI应该明显大于大气校正前的NDVI。分别在基于表观反射率的NDVI、基于本文模型得到的NDVI和基于USGS地表反射率产品的NDVI 3景影像上分别选取实际距离为5 km的截线，在截线上等距离取样本点，并进行统计计算分析，图5结果表明，基于USGS地表反射率产品的NDVI与基于本文模型的NDVI比较接近，两者的NDVI曲线明显高于表观反射率NDVI的曲线，这与上文的分析结果相一致。相关性分析表明基于模型的NDVI与基于USGS产品的NDVI的相关性（R²）最高达到了0.802 6。

阴影消除植被指数（SEVI）^［15］是一种在不借助DEM数据以及地面实测数据的条件下，能够有效地减弱山区中地形影响的植被指数，其采用红波段与近红外波段的非线性组合，增强山区中处于阴坡区域的植被指数值。SEVI计算如公式（3）所示。

其中：Bnir为近红外波段像元值，Br为红波段像元值，f∆为地形调节因子。

在SEVI的计算中，关键是f∆值的确定，该值的确定通过“寻优匹配”算法来实现，具体计算参见文献［16-18］。

由图7的7幅影像对比可以看到，经过本文模型校正前的NDVI影像以及RVI（比值植被指数）影像具有较为明显的地形纹理，处于阴坡处的区域相比于阳坡处的亮度较暗，这主要是由于地形影响，导致处于山体阴坡处的植被只能接收到较弱的太阳辐射而具有较低的亮度值，处于山体阳坡处的植被能接收到较强的太阳辐射而具有较高的亮度值。本文模型校正后的NDVI影像以及RVI影像的地形纹理有所消除，处于阴坡区域的植被信息得到了有效地增强，与表观反射率计算的SEVI在去除地形影响方面效果相当，USGS地表反射率产品也能够抑制地形影响，可能是本文模型和USGS的地表反射率产品算法均考虑了地形DEM的影响，在主要去除大气影响的同时具有一定的地形校正效果。

变量之间相关程度的指标由相关系数决定。相关系数的绝对值越大，变量的线性相关程度越高；绝对值越小，变量之间的线性相关程度越低。

评价地形校正效果的一个常用的客观指标是分析各植被指数与太阳入射角余弦值cosi的相关性，经过地形校正后，两者的相关性越低则校正效果越好。

表2为各植被指数与cosi的线性分析结果，如表2所示，本文模型校正后的NDVI-cosi线性拟合方程的斜率为-0.008 325 88，相关系数约为-0.010 168，SEVI-cosi线性拟合方程的斜率为-0.000 687 643，相关系数约为-0.005 353 70，校正前的NDVI-cosi线性拟合方程斜率为0.184 343，相关系数达到了0.256 224。相较于校正前NDVI-cosi线性拟合的结果，本文方法校正结果的斜率与相关系数都大幅降低，表明本文方法和SEVI的抗地形效果都较为明显，而USGS地表反射率产品在抗地形影响方面存在过校正现象（图8（f）~（g））。

（1）本文提出了一种基于查找表的Landsat-8 OLI遥感影像的大气校正方法，该方法由6S辐射传输模型生成查找表，输入的参数包括太阳天顶角、高程、水汽含量、臭氧浓度以及由MODIS数据反演得到的气溶胶光学厚度（AOT）。该查找表适用于所有的Landsat 8 OLI数据的地表反射率反演。该方法的创新之处在于采用了具有较高通用性的五维查找表，能够批量快速地生成地表反射率产品。

（2）通过光谱分析和交叉验证评价了模型计算结果的质量。通过目视分析，校正后的影像比原始影像更清晰。为了更准确地对模型质量进行评价，本文分别在基于表观反射率的NDVI、基于6S辐射传输模型的NDVI以及基于USGS地表反射率产品的NDVI 3景影像上分别选取实际距离为5 km的截线，在截线上等距离取样本点，并进行统计计算分析，通过分析发现基于6S辐射传输模型计算的NDVI比基于表观反射率的NDVI曲线更为符合典型植被的光谱特征。

（3）本文对不同算法减弱地形影响的效果进行了对比分析，通过对比分析基于表观反射率的NDVI、RVI、SEVI影像、基于6S模型的NDVI、RVI影像以及基于USGS地表反射率的NDVI、RVI影像，对-cosi植被指数进行线性相关分析，结果表明通过本文模型校正后的影像与SEVI在减弱地形影响方面的效果相近。