2013年4月26日发射的高分一号卫星(GF-1)，是中国重大科技专项“高分辨率对地观测系统”的首颗卫星。GF-1 PMS是搭载在GF-1卫星上的高空间分辨率的多光谱相机，该传感器具有高空间分辨率与高时间分辨率的特点，可以在国土、环境、农业等领域发挥重要作用。

大气校正是卫星数据应用中一个关键的前期处理工作，因为传感器接收到的辐射能量既包含地表信息也包含大气信息。大气的影响导致影像对比度降低，从而影响地表分类及定量信息的提取。

大气对遥感影像的作用主要包括分子散射，气溶胶散射以及水汽、臭氧、氧气、二氧化碳、气溶胶等的吸收作用(Liang等，1997)。对GF-1 PMS传感器而言，最主要的大气影响是由于气体分子引起的瑞利散射以及空气中气溶胶引起的米氏散射。瑞利散射影响相对较小，且可以依据经验高精度的计算，而气溶胶散射的影响则相对复杂，需要精确地提取气溶胶的光学厚度。因此，气溶胶光学厚度的精确提取是大气校正过程中首先要解决的一个关键问题。当前利用遥感影像提取陆地气溶胶光学厚度的算法已经发展了很多，而精度最高、可靠性最强的方式是浓密植被法(Kaufman和Sendra，1988，Kaufman等，1997a，Remer等，2005)。浓密植被法是利用陆地上浓密植被像元在红、蓝波段反射率较低，卫星传感器获取的辐射亮度值对气溶胶变化具有明显的指示作用的特点实现气溶胶的卫星遥感反演。浓密植被法气溶胶遥感的关键是地表反射率的精确确定。Kaufman和Sendra在1988年首先将浓密植被法应用于L and sat3 MSS传感器的气溶胶卫星遥感反演，经过近三十年的发展，浓密植被法地表反射率的确定经历了以下发展阶段：(1)基于植被指数的方法，其做法是根据NDVI确定浓密植被像元，对确定的浓密植被像元依据经验赋予地表反射率(Kaufman和Sendra，1988)，而由于NDVI随大气气溶胶的升高而显著降低(Kaufman等，1992)，致使在气溶胶光学厚度相对较大时，无法通过经验的阈值找到浓密植被像元，使得气溶胶反演无法实现，此外，不同的植被类型，其红、蓝波段的地表反射率有明显差异，依据经验对所有植被确定的统一的地表反射率精度较低；(2)基于2.1 μm波段的方法，Kaufman通过对地表测量的植被光谱以及航空数据获取的高分辨率植被像元的光谱特征的分析，发现植被在2.1 μm波段的反射率与红蓝波段的反射率具有稳定的相关关系(Kaufman等，1997a)，2.1 μm波段由于波长较长而受大气影响较小，利用2.1 μm波段受大气影响小且与红蓝波段存在线性关系的特点，确定浓密植被像元及其地表反射率，进而反演气溶胶光学厚度。暗目标法被广泛应用于MODIS(Kaufman等，1997b，Remer等，2005，Levy等，2007，Levy等，2013)、TM/ETM +(Liang等，1997，Ouaidrari和Vermote，1999)、VIIRS(Jackson等，2013)等含有短波红外波段传感器的气溶胶光学厚度的遥感反演和大气校正。

GF-1 PMS传感器缺少短波红外波段，导致大气气溶胶的高精度提取非常困难。针对缺少2.1 μm波段的传感器的气溶胶反演，孙林等(Sun等，2010，孙长奎等，2012)提出基于地表反射率数据支持的气溶胶反演方法，并利用MODIS的近红外波段支持HJ-1 CCD数据的气溶胶光学厚度反演和大气校正，取得了较好的效果。该方法通过对MODIS数据和HJ-1 CCD数据匹配后，为后者提供气溶胶遥感需要的地表反射率数据，然而异源数据精确配准是一项具有挑战性的任务，配准误差可能导致气溶胶的反演精度降低，甚至反演失败，限制了该算法的使用。为降低两类数据配准误差对地表反射率确定精度的影响，本文提出了使用区域NDVI分布百分比匹配的方法，回避了像元的直接匹配，将该方法用于对GF-1 PMS数据的气溶胶遥感和大气校正，有效提高了气溶胶反演的精度和大气校正的精度。

卫星传感器在可见光-近红外波段接收到的辐射能量可以表示为如下形式，

$\begin{align} & {{\rho }^{\text{*}}}\left(\tau，{{\theta }_{\text{s}}}，{{\theta }_{\text{v}}}，\varphi \right)={{\rho }^{a}}(\tau，{{\theta }_{\text{s}}}，{{\theta }_{\text{v}}}，\varphi)+\frac{{{\rho }_{\text{t}}}}{1-{{\rho }_{\text{t}}}S} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ T(\tau，{{\theta }_{\text{s}}})T(\tau，{{\theta }_{\text{v}}})\\ \end{align}$      (1)

式中，${\rho ^{\rm{*}}}\left({\tau，{\theta _{\rm s}}，{\theta _{\rm v}}，\varphi } \right)$为表观反射率，${\rho ^a}(\tau，{\theta _{\rm s}}，{\theta _{\rm v}}，\varphi)$ 为大气固有反射率(程辐射项)，${{{\rho _{\rm t}}}}$为目标反射率，S为大气半球反照率，$T(\tau，{\theta _{\rm s}})$、$T(\tau，{\theta _{\rm v}})$分别为大气下行透过率和大气上行透过率，${\theta _{\rm s}}$、${\theta _{\rm v}}$、$\varphi$分别为太阳天顶角、观测天顶角和相对方位角，τ为550 nm处气溶胶光学厚度。

气溶胶光学厚度反演和大气校正在本质上都是地气解耦合的过程。在式(1)中，除了气溶胶光学厚度和地表反射率之外，其他参数都可以通过卫星观测几何以及辐射传输方程中的大气模式、气溶胶模式等确定，即便如此，方程的个数仍然小于未知数的个数，为病态方程，要解算方程必须加入先验知识作为限制条件。目前气溶胶反演与大气校正的一般思路为：先以某些特殊地表类型的先验知识为限制条件，反演其上空气溶胶光学厚度，然后将反演的离散的气溶胶光学厚度插值、平滑用于大气校正。

浓密植被法气溶胶光学厚度反演，利用2.1 μm 短波红外波段受大气影响小且在浓密植被区域与红蓝波段地表反射率之间有稳定的比例关系的特点，来确定浓密植被区域红蓝波段的地表反射率，为式(1)提供了地表反射率数据，实现了浓密植被区域气溶胶光学厚度的反演。GF-1 PMS传感器缺少2.1 μm波段，但是根据浓密植被算法可以发现，浓密植被区域红光和蓝光波段的地表反射率存在一定的线性关系，

${\rho _{\rm red}} = k \cdot {\rho _{\rm blue}} + a$ (2)

式中，${\rho _{\rm red}}$和${\rho _{\rm blue}}$分别为红光和蓝光波段在浓密植被区域的地表反射率，k为红、蓝波段地表反射率比率，a为截距。

基于上述原理，首先确定GF-1 PMS数据浓密植被像元的分布，然后根据GF-1 PMS数据浓密植被像元红、蓝波段地表反射率之间的线性关系，结合式(1)和式(2)，解算辐射传输方程即可获得浓密植被像元气溶胶光学厚度。

在缺少短波红外波段的条件下，寻求一种稳定可靠的方式来确定浓密植被像元非常困难。Kaufman等利用原始影像的NDVI以及近红外波段的表观反射率来确定浓密植被像元(Kaufman和Sendra，1988)，然而这些参数本身就受气溶胶光学厚度等大气条件的影响，难以通过确定一个固定阈值，来准确确定影像中的浓密植被像元。

本文提出使用现有的地表反射率数据集辅助提取GF-1 PMS的浓密植被像元，使用的数据集为MODIS的地表反射率产品数据。为降低空间匹配对浓密植被像元提取产生的误差，这里提出了利用浓密植被“百分比匹配”的方式获得GF-1PMS浓密植被像元的分布。所谓百分比匹配是指在MODIS和GF-1PMS数据区域相同的前提下，首先利用MODIS地表反射率数据统计出该区域浓密植被百分比P_DDV；然后将GF-1PMS传感器多光谱数据重采样为与地表反射率数据集相同的空间分辨率，利用重采样后的表观反射率影像计算NDVI；最后根据P_DDV确定PMS数据中浓密植被像元NDVI 的阈值，大于该阈值的确定为浓密植被像元，参与之后的气溶胶光学厚度反演。

对于红、蓝波段地表反射率之间关系的确定，是利用LOPEX’93以及USGS波谱库中的植被光谱曲线与GF-1 PMS传感器多光谱波段的光谱响应函数模拟计算，获得各种植被红、蓝波段的地表反射率，并对模拟的地表反射率进行统计分析。图 1为回归分析结果，GF-1 PMS传感器在浓密植被区域红、蓝波段地表反射率的关系为：${\rho _{\rm red}} $= 1.2301+${\rho _{\rm blue}}$+0.0022。

图 1 GF-1 PMS传感器红光波段与蓝光波段线性回归分析

Figure 1 Linear regression analysis of GF-1 PMS band 1 and band 3

基于确定的浓密植被像元和拟合的红蓝波段的相关关系，可以实现GF-1 PMS数据的气溶胶光学厚度遥感反演。反演过程中辐射传输方程的求解涉及的参数较多，过程复杂，这里使用查找表的方法求解，利用6S辐射传输模式(Second Simulation of Satellite Sensor in Solar Spectrum)(Vermote等，1997)构建查找表，构建查找表涉及的参数包括成像几何、气溶胶光学厚度、程辐射项、大气半球反照率以及大气透过率等。

GF-1卫星上搭载了两个PMS相机(PMS1、PMS2)，相机获取的数据包括2 m全色图像和8m多光谱图像，其参数设置如表 1所示。两台PMS相机组合使传感器在空间分辨率高的同时具备较大的幅宽，另外传感器具有良好的侧摆成像能力，可以在应急状态下对特定的区域进行高时间频次观测。GF-1卫星轨道高度约为645 km，PMS相机单景幅宽约30 km，通常情况下其观测天顶角小于10°，接近垂直观测。

表 1 GF-1 PMS参数设置
Table 1 The Parameters settings of GF-1 PMS

下载CSV 

参  数	PMS相机
光谱范围	全色	0.45—0.90 μm
	多光谱	0.45—0.52 μm
		0.52—0.59 μm
		0.63—0.69 μm
		0.77—0.89 μm
空间分辨率	全色	2 m
	多光谱	8 m
幅宽	60 km(2台相机组合)
重访周期(侧摆时)	4 d
覆盖周期(不侧摆)	41 d

PMS相机多光谱数据有4个波段，为常见的蓝(band 1，0.45—0.52 μm)、绿(band 2，0.52—0.59 μm)、红(band 3，0.63—0.69 μm)、近红外(band 4，0.77—0.89 μm)4个可见光/近红外波段。两个相机的光谱响应函数几乎相同，其多光谱波段的光谱响应函数如图 2所示。

图 2 GF-1 PMS相机光谱响应函数

Figure 2 The spectral response function of GF-1 PMS

本文利用MODIS/Terra 8天合成全球地表反射率产品MOD09A1(正弦投影，空间分辨率500m，同一轨道号每年有46景影像)构建地表反射率数据集。

从美国国家航空航天局NASA地球观测系统的数据和信息系统(EOSDIS)下载2008年—2012年5年间覆盖中国大部分陆地区域的所有8天合成的MOD09A1地表反射率数据。对下载的MOD09A1数据进行投影转换、镶嵌等预处理，再利用全国矢量图进行裁剪获得5年全国陆地地表反射率数据集。为了降低云覆盖的影响，假设5年内地表类型变化不大，按最小值合成的方式将5年间同一日期的5景影像合成为一景影像，该影像基本全部由晴空像元组成，最终获得46景覆盖中国大部分陆地区域的地表反射率数据集。图 3是第225—232天的5年合成地表反射率数据集彩色合成影像，从影像中可以看出，经过最小值合成后有效消除了云像元的影响，整景影像基本全部由晴空像元组成，可以通过该数据集获得指定时间、空间的地表反射率信息。在该地表反射率数据集的支持下，可以利用“百分比匹配”的方式获得GF-1 PMS传感器多光谱影像中浓密植被像元的分布。

图 3 第225—232天5年合成地表反射率数据彩色合成影像

Figure 3 The false color image of the 5-year composite surface reflectance for the 225th to 232th days each year

为了评价气溶胶反演算法的精度，获取了北京(116.381°E，39.977°N)和太湖(120.215°E，31.421°N)两个全球气溶胶自动观测网(AERONET)站点的地基观测数据进行对比验证。AERONET(Holben等，1998，Holben等，2001)是由NASA和法国国家科学研究中心CNRS共同建立的一个全球范围的地基气溶胶观测网络，目前共有大约800个观测站点(陈好等，2013)。站点每隔15 min获取一次数据，可以同时获得7个波段(0.340—1.020 μm)的气溶胶光学厚度值，精度大约为0.01—0.02(Holben等，2001)，这一精度是卫星反演精度的3—5倍(Remer等，2009)，Level2.0产品经过云滤除和质量控制(Smirnov等，2000)，常被用于卫星观测气溶胶光学厚度的精度验证。观测网络所有站点数据向全球用户免费共享。

为了保证卫星遥感反演值与站点观测值在光谱、时间、空间上的一致性，需要做出一些合理的假设，并对两种数据源进行适当的处理。在光谱尺度上，由于卫星遥感反演结果为550 nm处气溶胶光学厚度，而AERONET站点没有观测这一波长的气溶胶光学厚度，因此需要基于Angstrom指数关系(Ångström，1929)，利用AERONET数据中440 nm与870 nm波长对应的气溶胶光学厚度计算550 nm波长处的气溶胶光学厚度。

$\tau \lambda \text{=}\beta {{\lambda }^{-\alpha }} $ (3)

式中，λ为波长；β为大气浑浊度系数，与气溶胶粒子总数、粒子谱分布和折射指数有关；α为Angstrom指数，反映粒子大小，值越大气溶胶粒子越小。

在时间尺度上，认为气溶胶在同一地区短时间内是线性变化的，基于这一假设选择与遥感影像成像时刻最接近的前后两个时刻的站点观测值，对它们线性插值作为观测站点在遥感影像成像时刻的站点观测值。

$\tau \text{=}{{\tau }_{1}}\frac{{{t}_{2}}-t}{{{t}_{2}}-{{t}_{1}}}+{{\tau }_{2}}\frac{t-{{t}_{1}}}{{{t}_{2}}-{{t}_{1}}} $ (4)

式中，τ为卫星图像成像时刻的气溶胶光学厚度，t为成像时刻，t₁、τ₁ 为成像前站点观测值的观测时间与气溶胶光学厚度值，t₂、τ₂为成像后站点观测值的观测时间与气溶胶光学厚度值。

根据两类数据的空间差异以及气溶胶空间分布的特点，将GF-1 PMS数据转换到500 m分辨率的尺度进行气溶胶光学厚度的提取。本文选用了5景GF-1 PMS多光谱遥感影像进行了气溶胶遥感反演实验，影像的缩略图和反演的气溶胶光学厚度图像如图 4所示，其中黄色十字星代表AERONET站点的大致位置。反演结果与AERONET站点气溶胶光学厚度Level 2.0产品以及MODIS气溶胶产品MOD04进行对比，表 2为对比结果。对比时，假定气溶胶光学厚度在较小的空间范围内变化平缓(Liang等，2001)，选取与AERONET站点对应的遥感影像上21 × 21个像元反演均值作为GF-1 PMS反演结果与其对比。

图 4 实验所用的5景影像的缩略图以及AOD反演结果

Figure 4 Thumbnail images and AOD results of the five scenes used in the experiments

表 2 GF-1卫星AOD反演值与AERONET站点观测值对比结果
Table 2 Results of GF-1 PMS AOD against AERONET AOD measurements

下载CSV 

站点	成像日期	AERONET值	MOD04结果	反演结果	绝对误差	相对误差
北京站	2013-06-19	0.18247	0.4840	0.26647	0.08400	0.46038
北京站	2013-06-27	0.32525	0.8610	0.54453	0.21929	0.67422
北京站	2013-08-10	0.95988	1.3880	0.94948	–0.01041	0.01084
北京站	2013-08-14	1.20639	1.2660	1.06020	–0.14623	0.12121
太湖站	2013-07-24	0.39612	0.3805	0.40000	0.00388	0.00980

虽然受数据源的限制只获得了5景影像实验，但通过分析这5景影像的反演精度，可以发现几个非常明显的规律：首先，气溶胶光学厚度较高时(图 4(c)(d))，卫星反演值低于AERONET站点地面实测值，气溶胶光学厚度较低时(图 4(a)(b))，卫星反演值高于AERONET站点地面实测值，但总体上气溶胶光学厚度较低时反演的相对误差较大。这一方面是由于AOD 较小的情况下，较小的绝对误差也会产生较大的相对误差，另一方面也可能是由于气溶胶光学厚度较小的情况下，卫星接收的信号中气溶胶的贡献占得比例较小，导致卫星图像对气溶胶光学厚度的指示作用降低(Bréon等，2011，方莉等，2013)。其次，云影响气溶胶光学厚度的反演精度，尤其是厚云影响明显。从5景图像中可以看出，图 4(b)(d)影像中云层厚度较大，且遮盖于浓密植被上方，影响了DDV像元的确定，从而影响最终AOD的反演精度。而图 4(c)影像虽然有云的影响，但云层相对较薄，不影响DDV像元的确定，最终AOD的反演精度较高。再次，地形可能是影响AOD反演精度的重要因素。对于图 4(a)(b)影像，AOD反演精度较低的原因除了AOD较小外，也有地形的影响。图 4(a)(b)影像中，DDV像元主要分布在香山区域，而山区地形起伏较大，阴影、双向反射率特征都可能对地表的光谱特征产生影响。本文采用地表为朗伯体的假设反演AOD，在山区这一假设受到很大的限制，由此导致AOD反演的精度偏低。最后，地表反射率对反演AOD的影响也在5景影像中有所体现。北京地区AOD的反演精度低于太湖地区AOD的反演精度，这是由于北京地区地表反射率较高，大气对卫星图像的贡献相对较小，且高反射率的城市地区对浓密植被的邻近效应明显，导致AOD反演精度较低。对于图 4(e)影像太湖地区则是一个反演AOD的理想区域：该地区浓密植被丰富，适宜反演AOD的DDV像元多；水体和浓密植被是主要地表类型，地表覆盖较为均一，邻近效应小，地表反射率低，大气对卫星信号的贡献大；陆地以平原和低矮丘陵为主，地形变化平缓；另外，该景图像AOD大小适中，且不受云的影响。在以上因素的共同作用下，该景图像反演AOD的精度极高。

图 5为AOD的反演结果与AERONET站点观测值的回归分析结果。AOD卫星反演值与站点观测值的回归方程为y=0.8593x+0.0905，相关系数为0.92，均方根误差为0.115。Kaufman等人(1988)指出，气溶胶光学厚度误差与大气校正误差的对应关系约为Δτ~10Δρ，即由气溶胶光学厚度误差导致的大气校正误差大约是气溶胶光学厚度误差本身的1/10左右。从总体上看，5景影像中有4景AOD的反演精度优于0.15，由此导致的大气校正误差小于0.015，具有良好的效果。

图 5 GF-1 PMS气溶胶光学厚度反演值与AERONET站点观测值以及MOD04对比结果图

Figure 5 GF-1 PMS AOD against AERONET AOD measurements and MOD04

将反演的AOD结果影像插值、平滑，获得每个像素的AOD值，再将该值连同表观反射率、观测几何、大气模式等信息一起代入式(1)中即可计算每个像素的地表反射率值，从而实现大气校正。接下来分别从视觉效果、定性统计、定量分析的角度对大气校正的效果进行评价。

选取2013-08-14北京地区GF-1 PMS传感器多光谱遥感影像，利用本文提供的大气校正方法对影像进行大气校正， 图 6是实验结果的截图，影像采用标准假彩色显示。通过对比大气校正前后的影像可以看出，影像“模糊”程度明显降低，纹理细节更加清晰，影像对比度有显著提高。

图 6 大气校正前后影像按标准假彩色显示效果对比

Figure 6 The standard false color images before and after the atmospheric correction

对校正前后的影像截图部分进行基本统计和直方图统计，统计结果如表 3和图 7所示。从统计结果可以看出，大气校正后的影像4个波段的标准差都明显提高，这与目视效果完全吻合，说明影像的对比度有显著提高。另外，观察大气校正前后的直方图可以发现，大气校正后影像的直方图“平滑”了许多，说明大气校正过程中滤除了遥感影像中的大气噪声。

表 3 大气校正前后影像截图部分的均值和标准差
Table 3 Statistical values of the selected part before and after atmospheric correction

下载CSV 

波段	校正前均值	校正后均值	校正前标准差	校正后标准差
第1波段	0.1606	0.0637	0.0072	0.0203
第2波段	0.1337	0.0493	0.0087	0.0232
第3波段	0.1111	0.0567	0.0123	0.0271
第4波段	0.1437	0.1805	0.0235	0.0463

图 7 大气校正前后各波段灰度直方图对比

Figure 7 The histogram comparison of the four bands

选取2013-06-26敦煌辐射定标场的两景GF-1 PMS传感器多光谱影像。两景影像中一景分布有石膏矿，代表地表反射率高的地物，如图 8(a)所示；另一景中分布有水体，代表地表反射率低的地物，如图 8(b)所示。分别用本文算法对两景影像进行大气校正，将石膏矿和水体两种地物的表观反射率、地表反射率大气校正反演值、地表反射率地面同步实测值进行比较，比较结果如表 4和图 9所示。

图 8 目标点周围环境

Figure 8 Targets surroundings

图 9 石膏矿与水体大气校正前后反射率与地面实测地表反射率波谱曲线对比图

Figure 9 Comparison of the apparent, corrected and the field measured reflectance

表 4 石膏矿与水体大气校正前后反射率与地面实测地表反射率对比
Table 4 Results of the apparent, corrected and the field measured reflectance

下载CSV 

		第1波段	第2波段	第3波段	第4波段
	表观反射率	0.2623	0.2813	0.3154	0.3290
石膏矿——高亮地表	地表反射率反演值	0.2793	0.3236	0.3725	0.4028
	地面实测地表反射率	0.2716	0.3272	0.3755	0.4004
	绝对误差	0.0077	–0.0036	–0.0030	0.0024
	表观反射率	0.0932	0.0777	0.0466	0.0329
水体——暗地表	地表反射率反演值	0.0433	0.0487	0.0297	0.0270
	地面实测地表反射率	0.0345	0.0438	0.0243	0.0160
	绝对误差	0.0088	0.0049	0.0054	0.0110

从表 4中可以发现，水体在各个波段的地表反射率反演值都小于其表观反射率值；石膏矿则正好相反，在各个波段的地表反射率反演值都大于其表观反射率值。这与大气可以提高低地表反射率地物的表观反射率，降低高地表反射率地物的表观反射率的结论(Fraser和Kaufman，1985)是吻合的。观察图 9可以发现，无论从走势还是接近程度看，大气校正后地表反射率值较之表观反射率更接近于地表测量的反射率。实验结果表明，除了水体在近红外波段的地表反射率反演误差略大于0.01以外，其他地表反射率的反演误差都在±0.01以内，大气校正的精度基本可以达到Δρ≤0.015。