0 引言

遥感技术的应用越来越广泛，而随着遥感应用的定量化、精细化发展，越来越多的定量化遥感涉及到光学分析^{[1, 2, 3, 4, 5]}。因传感器获得的信号是地面信号和大气信号的叠加，大气影响不仅造成成像传感器图像模糊、对比度下降、信噪比降低以及细节损失^[6]，还使得遥感数据在后续的植被指数、地表温度、叶绿素浓度等定量化反演中所获得的结果严重偏离实际值^[7];因此需要重视大气对遥感数据的影响。

大气校正是遥感辐射校正的重要步骤，是由遥感影像获取反映地表真实反射率的一个必不可少的环节，对于定量遥感而言尤为重要^{[8, 9]}。这是由于传感器在获取信息过程中受到大气分子、气溶胶和云粒子等大气成分吸收与散射的影响，使其获取的遥感信息中带有一定非目标地物的成像信息，数据预处理的精度达不到定量分析的高度。一方面，大气的吸收、散射作用削弱了地表反射能量，对图像信号有衰减作用；另一方面，大气程辐射对图像信号有增加作用^[10]。消除这些大气影响的处理，称为大气校正^[11]。

以往进行大气校正过程中，假设大气条件均一，整幅影像仅使用一套参数进行校正。但在实际情况下，卫星观测地面生成的影像每一个像元的大气参数、观测角度、海拔高度等参数均不相同，特别是在雾霾天气的情况下，气溶胶分布的不均匀更较为直观地影响着卫星影像大气校正的结果。

因此综合上述情况，笔者尝试逐像元考虑各方面参数进行大气校正，选用MODIS(moderate-resolution imaging spectroradiometer)数据作为研究对象，基于6S辐射校正模型统计分析各个参数对反射率的敏感性，生成查找表，逐像元考虑大气、角度、高程等参数进行大气校正，以提高大气校正的精度。

6S模型由Eric Vermote等人在5S(the simulation of the satellite signal in the solar speetrum radiative code)模型的基础上进行改进，采用SOS(successive order of scattering)的散射计算方法，建立在辐射传输理论基础之上，准确模拟太阳——目标物——传感器路径上的大气影响，应用范围广，受研究区特点及目标类型等的影响较小，精度高，与LOWTRAN、MODTRAN等模型同为目前发展比较成熟的大气订正模型^[12]。6S模型适用于0.25~4.00 μm波长范围内电磁波的大气辐射传输的模拟。利用该模型进行大气校正的工作流程是^[13]：按顺序将几何参数、大气模式、气溶胶模式、可见度、海拔高度、观测波段等参数输入6S模型;通过计算模拟太阳辐射在大气中的辐射状况;给出大气校正系数X_a、X_b和X_c^[14]。

假设目标表面为均匀朗伯体表面，不考虑气体吸收，在6S模型中卫星观测的表观反射率可用式(1)^[15]表达：

点击浏览公式

式中：ρ_TOA为大气表观反射率；ρ₀为大气的路径辐射项等效反射率；ρ_S为地表反射率；θ_s为太阳天顶角；θ_v为卫星天顶角；φ_s为太阳方位角；φ_v为卫星方位角；φ_s——φ_v代表相对方位角；S代表大气球面反照率；T(θ_s)代表大气下行透过率；T(θ_v) 代表上行辐射总透过率。

运行6S辐射传输模型后，可得到大气校正系数X_a、X_b、X_c，然后按式(2)计算得到地表反射率^[16]：

点击浏览公式

式中:Y为中间变量;L为定标后的辐射亮度值。

本文选用的MODIS遥感影像以安徽省为例，区域范围为东经114.9°——119.6°，北纬29.4°——34.6°，成像日期为2013年7月3日。

根据研究区遥感影像图，6S模型的初始输入参数见表 1，同时假定地表具有均一的朗伯反射特性。其中，气溶胶厚度利用MOD04气溶胶产品得出，地面高程利用SRTM(shuttle radar topography mission)地形产品数据得出，二者均经过投影、剪裁、重采样、辐射定标等预处理。

对于处理大范围遥感影像图，参数一致的假设显然是导致地表反射率校正误差的原因^[17]。因此通过对敏感性的分析，我们能更深层次地了解各种参数对大气校正结果的影响程度，为大气校正参数选取提供参考。根据6S模型模拟，影响大气校正获得地表反射率的参数包括几何参数、目标和传感器的高程参数、大气模式、气溶胶模式以及光谱条件。因此，笔者选择了太阳天顶角、传感器天顶角、相对方位角、大气气溶胶厚度、观测波段以及地表海拔这6个参数的变化来分析不同参数对大气校正的影响。利用IDL语言线下调用6S辐射模型生成查找表，统计分析各参数关于地表反射率敏感度的优劣(图 1)，各参数的取值范围如下：太阳天顶角0°~70°，卫星天顶角0°~80°，气溶胶光学厚度0.00~1.95，海拔高程0.0~1.9 km。

图 1 不同观测波段上的气溶胶光学厚度(a)、太阳天顶角(b)、卫星天顶角(c)、海拔高程(d)及相对方位角(e)对反射率的敏感度 Fig.1 AOD(a), solar zenith(b), sensor zenith(c), ground altitude(d) and relative azimuth(e) of the sensitivity to the reflectivity on different observation bands

图选项

据图 1分析，气溶胶光学厚度对反射率的影响稳定且明显，气溶胶厚度越大对反射率的敏感度也就越好(图 1a)。太阳天顶角在角度小于25°时，对地表反射率的影响不明显，而太阳天顶角大于25°时，太阳天顶角对反射率的敏感度明显增高，甚至高于气溶胶光学厚度对于反射率的敏感度(图 1b)。卫星天顶角在MODIS的波段3和波段4的敏感度明显高于其他波段，随着角度的增大反射率略有下降，卫星天顶角整体敏感度不高(图 1c)。海拔高度对反射率的敏感度同样较低，笔者仅计算了海拔高度在1.9 km以下的地表高程，海拔高度仅在MODIS的波段3和波段7对反射率有一定的影响，其他波段的敏感度几乎可忽略不计(图 1d)。另外，笔者也注意到当太阳天顶角变化较大时与变化较小时观测的海拔高程对反射率的敏感度几乎一致。而相对方位角对反射率的敏感度最差，相对方位角对应的大气校正系数仅X_b发生细微变化(图 1e)；根据公式(2)、(3)可得出相对方位角对反射率的影响几乎可忽略。

对地表反射率影响敏感度较高的有气溶胶光学厚度及太阳天顶角，而相对方位角敏感度最差。经过针对上述各个参数的敏感性分析，我们在对某一景影像做大气校正时，可以根据这景卫星图像的各参数情况，选择性地忽视一些敏感度较低的参数，这样在提高大气校正效率的同时也一定程度地缩小了大气校正的工作量，同时也扩展了利用辐射传输模型进行大气校正的应用。

笔者在对研究区做逐像元大气校正时，考虑到因参数较多导致数据量偏大、耗时较多，因此针对研究区利用IDL语言调用6S辐射模型做查找表，利用6S模型离线计算不同太阳天顶角、卫星天顶角、气溶胶光学厚度以及地表海拔高程情况下的校正系数X_a、X_b、X_c，生成研究区上述各个参数范围内的大气校正系数。由于大气校正是逐像元分析而且又加入考虑了多个影响大气反射率的参数，所以针对工作区做查找表的好处在于，输入大气校正参数的范围被缩小，根据各参数范围分别调整步长，这样在缩小工作量的同时可以更加精细化地计算大气校正的参数。

本文研究区面积较小，研究区内参数的范围比较狭小，所以本文参数的取值范围如下：太阳天顶角为0°~17°时步长设为3°，太阳天顶角为24°~60°时步长设为12°，共11个值；卫星天顶角0°~75°以15°为步长，共6个值；气溶胶光学厚度0~1.95，步长不规则，共11个值；地表海拔0~1.9 km，步长不规则，共9个值；以及MODIS的前7个波段。将上述参数组合代入6S辐射模型，得到不同参数下大气校正系数的查找表。

逐像元大气校正步骤：

1)根据遥感图像的成像时间，获取研究区那一时刻的大气参数、角度信息以及高程数据，分别进行投影、剪裁、重采样等预处理，使得获取的各类参数对应于遥感影像信息。

2)利用获取的研究区内的各类参数信息，分别取均值作为输入6S辐射模型的初始参数，得到6S辐射校正输出表。

3)在辐射校正输出表中确定需要参数的位置，根据研究区各参数范围分别设计步长，并利用IDL编程循环调用6S辐射校正模型，以获取不同参数不同组合的大气校正系数查找表^[14]。

4)对大气校正系数查找表进行线性插值，并根据研究区各像元的大气、角度、高程等信息查询查找表中对应的大气校正系数，根据像元对应的大气校正系数计算该像元处的地表反射率。

笔者将查找表逐像元大气校正所得出的结果与统一参数大气校正以及校正前的表观反射率结果进行了对比，采用人工目视解译的方法随机选取了陆地和水体做为典型地类分别获取像元处的地表反射率，结果如图 2所示。由于气溶胶及观测角度等参数的不同，统一参数大气校正与逐像元大气校正的精度差值比较明显。

a.陆地；b.水体；c.频率。 图 2 三种不同方法的校正结果在各观测波段上的比较 Fig.2 Comparison of the correction result of three method for each observation band

图选项

校正前卫星定标后得到的表观反射率虽然消除了太阳天顶角、日地距离和太阳辐射量的差异，但是不能消除由大气程辐射、瑞利散射和吸收带来的影响^[18]。这些影响造成的大气衰减使得卫星传感器接收到的信号很大部分不是地面反射率的光谱信息，校正前后差值很大。地表反射率图像中地物在可见光波段(波段1——4)反射率减小，其中蓝光波段(波段3)的瑞利散射最强，反射率减小最为明显，绿光(波段4)次之(图 2a、b)。

根据研究区校正前后反射率频率变化图(图 2c)可以看出,校正后的反射率峰值向低反射率方向移动，整体上校正前的表观反射率值大于校正后的反射率值。

由于6S模型计算耗时太长，取研究区像元第497行来进行验证。图 3给出了研究区第497行扫描线上自西向东各个像元的气溶胶光学厚度、太阳天顶角、卫星天顶角、海拔高程和相对方位角的情况，结合各参数对反射率的敏感度，对比校正前后反射率的变化情况。本文以红光波段为例，根据扫描线上的像元参数进行逐像元大气校正，将校正结果与统一参数大气校正以及表观反射率的结果进行对比(图 4a)。

图 3 扫描线上各像元的气溶胶光学厚度(a)、太阳天顶角(b)、卫星天顶角(c)、海拔高程(d)、相对方位角(e)情况 Fig.3 AOD(a), solar zenith(b), sensor zenith(c), ground altitude(d) and relative azimuth(e) of the pixels on the scan line

图选项

大气校正削弱了大气瑞利散射和气溶胶散射所导致的红光波段信号增强，而地形因素又使得校正后的反射率大于表观反射率，校正后的反射率与实际情况更为接近。利用MODIS地表反射率产品对大气校正的结果进行验证，统一参数、逐像元大气校正的相对误差分别控制在26.9%和12.7%以内。因此，对于影像范围较大，或气溶胶和地面高程分布不均匀的山区或者农田等地逐像元大气校正方法更具优势。

以植被指数为例证明大气校正精度对后续遥感定量计算的影响程度^[19](图 4b)。从图 4b可以看出，校正前后同一剖面(第497行)两条NDVI曲线的分布和变化趋势吻合，峰谷和峰值基本一致，逐像元比单一参数大气校正方法计算的NDVI平均高出14.4%。究其原因，经过大气校正后，可见光的反射率降低，近红外波段的反射率升高，两波段的对比度增高，光在大气传播过程中所受到的衰减得到弥补，则影像的植被指数呈现增加的趋势^[18]。

图 4 第497行反射率(a)、植被指数(b)在三种大气校正方法下的结果对比 Fig.4 The 497 line reflectivity and vegetation index results were compared in three atmospheric correction methods

图选项

分析了6S模型中各种参数对反射率的敏感度，并详细叙述了逐像元大气校正方法。结论如下：

1)利用查找表逐像元大气纠正的算法，根据各个像元的大气状况、观测角度、地面高程等条件差异，综合考虑了各种敏感度不同的参数，与以往常用的基于整景遥感影像获取平均参数并进行大气校正的算法相比，逐像元大气校正算法与实际情况较为接近，能够更加精确地对遥感影像进行大气校正并获取地物的真实反射率。

2)利用MODIS地表反射率产品对大气校正的结果进行验证，单一参数、逐像元大气校正的相对误差分别控制在26.9%和12.7%以内；在以植被指数(NDVI)为例的后续遥感定量化反演过程中，逐像元比单一参数大气校正方法计算的NDVI平均高出14.4%。

3)查找表的建立是大气校正精度提高的关键，在一定范围内，参数的步长过长影响校正的精度，而步长过小则使得数据量过大，冗余数据较多。因此可以根据研究区情况设定参数的范围，并分别设计各参数的步长，根据研究区来建立查找表，而不是固定用一套查找表。

4)经过大气校正后的植被指数整体上明显提高，本文的逐像元大气校正方法可以用于分析影像范围较大，或气溶胶和地面高程分布不均匀的山区或者农田等地。