资源三号卫星(以下简称ZY-3卫星)是我国第一颗民用高分辨率测绘卫星，其多光谱数据空间分辨率为5.8 m，量化级别为16 bit，包括3个可见光波段(0.45~0.52、0.52~0.59、0.63~0.69 μm)和1个近红外波段(0.77~0.89 μm).ZY-3多光谱图像在分辨率、量化级别、清晰度、色调及纹理特征方面都优于SPOT5图像，能够在对地观测的各个领域发挥重要的作用(李霖等，2014).然而，目前关于ZY-3卫星的研究大多集中于测绘行业，对于其他行业定量遥感的研究很少.传感器接收到的信息是大气与地表信息的耦合，提取气溶胶信息进行大气纠正是定量遥感必须考虑的因素，也是气溶胶遥感反演的核心.因此，针对ZY-3数据特点开发提取气溶胶信息的算法具有以下3方面意义：①ZY-3高空间分辨率特性更适合于城市研究尺度，对于城市气溶胶的空间分布、污染源及区域输送具有重要意义；②ZY-3能够为缺少光学数据的多云多雨地区提供新的气溶胶遥感数据源；③从其他卫星获得的气溶胶数据分辨率偏低，不能很好地为ZY-3大气校正服务，利用ZY-3数据本身获得气溶胶进行大气纠正是数据定量化应用的基础.

目前，监测陆地气溶胶的方法主要有暗像元法(Kaufman et al., 1988)、结构函数法(Tanré et al., 1988)、基于地表反射率数据库深蓝算法(Hsu et al., 2004)、双星协同法(Tang et al., 2005)、BAER算法(von Hoyningen-Huene et al., 2003; 2006)等.除暗像元法及BAER算法外，其余算法需要额外的卫星观测数据作为辅助实现地气解耦.BAER算法依靠区域长期大气观测数据建立气溶胶模式，难以在缺乏地基观测数据的地区使用.考虑到气溶胶在时空上变化万千，最可行的方法就是从图像自身获取成像当时的气溶胶参数，同时解算大气气溶胶和地表反射率信息(Liang et al., 1997).

浓密植被暗像元法(DDV)就是一种经典且常用的基于图像自身的气溶胶提取算法.然而，ZY-3卫星缺少短波红外波段，不适用基于MODIS数据提出的DDV法.杨磊等(2013)通过结合MODIS数据反演的气溶胶光学厚度(1 km)，基于6S辐射传输模型构建查找表，对ZY-3卫星多光谱数据进行了大气校正.实际上，目前已有许多利用浓密植被红蓝波段进行气溶胶反演的算法.Richter等(2006)修正了浓密植被暗像元的判定条件及ETM红波段与近红外波段地表反射率线性关系，使其能够适用于只含蓝-绿-红及近红外4个波段的卫星数据.何颖清等(2010)提出了一种基于红蓝波段自动选取多个暗 像元的 TM 影像大气纠正方法，但模型假设气溶胶 呈瑞利散射建立，与实际存在一定误差.Gilabert等(1994)提出了一种仅利用程辐射信息估算气溶胶光学参数并直接计算AOD的方法，因方法适用性较强，许多学者在此基础上进行了改进.Zhang等(2005)在Gilabert模型中增加了气溶胶垂向分布的考虑，但其假设区域大气状况相同，无法反映大气状况复杂的城区AOD分布.宋巍巍等(2008)利用Gilabert模型，选择多个暗像元，通过空间插值的方法获得广州市气溶胶AOD分布.Themistocleous等(2013)将Gilabert模型中的指数项简化为二次项，反演获得利马索尔城区AOD分布.然而，Gilabert模型未考虑暗像元地表漫反射信息，与实际具有一定的误差.本文在此方面进行改进，旨在开发一种仅利用图像本身信息自动反演AOD分布，而精度尽可能高的算法.由于目前关于ZY-3之类高分辨率数据气溶胶反演及大气纠正的研究较少，本文将算法用于分析ZY-3数据香港地区气溶胶AOD分布状况，识别局地大气污染源，以期为高空间分辨率数据应用于城市区域气溶胶自动监测提供参考.

香港及其所处的珠三角地区的大气污染问题一直受到密切关注，而气溶胶信息在城市亮地表区域较难剥离，因此，对香港这类大都市的气溶胶遥感研究具有现实和科学意义.另一方面，香港地区虽然范围较小，但其地形起伏较大，气溶胶在三维空间内变化较大，适宜于检验算法的有效性.本文采用的香港地区ZY-3多光谱影像从中国资源卫星应用中心网站获得，其成像日期为2013年3月8日.辅助数据为1 ∶ 5万DEM，在http://www.gscloud.cn/网站下载.用于验证的MODIS气溶胶C051产品数据在http://ladsweb.nascom.nasa.gov网站下载获得.

ZY-3遥感影像原始观测数据为DN值，根据公式(1)转化为行星反射率ρ_TOA(λ)：

式中，θ_s为太阳天顶角(°)，d是日地距离修正因子，可查表获得，Gain为增益，E₀为大气层外太阳辐照度，取值如表 1所示.

表1(Table 1)

表1 ZY-3多光谱数据辐射纠正参数及大气层外太阳辐照度 Table 1 ZY-3 MUX calibration coefficients and solar irradiance in the upper atmosphere

表1 ZY-3多光谱数据辐射纠正参数及大气层外太阳辐照度 Table 1 ZY-3 MUX calibration coefficients and solar irradiance in the upper atmosphere 波段 Gain /(W·m－2·sr－1·μm－1) E0/(W·m－2·μm－1) B and 1 0.2551 2012 B and 2 0.2353 1850 B and 3 0.1944 1550 B and 4 0.2107 1100

假设天空辐照度各向同性且地面为朗伯面反射，忽略地表环境光，仅考虑一次散射，则卫星观测的行星反射率ρ_TOA(λ)可以近似地表示为程辐射反射率及地表漫反射两部分.通常程辐射可分离为大气分子瑞利散射反射率ρ_pr(λ)和气溶胶米散射反射率ρ_pa(λ)两部分(Gordon，1978).地表漫反射由辐射来源可分为太阳直射光反射率ρ_b(λ)、大气分子下行散射光反射率ρ_dr(λ)及气溶胶下行散射光反射率ρ_da(λ)3部分(Iqbal，1983).则行星反射率ρ_TOA(λ)可表示为以下5部分组成：

ρ_pr(λ)的表达式为(Saunders，1990)：

大气分子瑞利散射各向同性，其散射光下行分量的比例为1/2.对于非各向同性的气溶胶散射，可用F_c(θ_s)来描述散射光下行部分与总能量的比例.Iqbal(1983)给出了不同太阳天顶角下F_c(θ_s)的离散值，据此通过线性插值计算出成像时刻的F_c(θ_s)值.

由于辐射传输方程的形式复杂，直接从数值上分离大气和地表信息是难以实现的.可以寻找影像上一些地表反射率很低的暗像元，先近似认为ρ_s0，忽略地表漫反射，利用程辐射信息估算气溶胶散射光学参数P_a和ω₀.进而采用迭代算法考虑地表漫反射，求取更接近实际的τ_a值.最后将这些暗像元的大气光学参数进行插值，从而实现整景图像的地气解耦.

浓密植被在红、蓝波段的地表反射率很小，是一种典型的暗像元.而阴影植被的太阳直射光几乎为零，较浓密植被更适宜作为暗像元，可认为其行星反射率与程辐射反射率近似相等.假设气溶胶呈Junge分布，根据ngström公式气溶胶光学厚度可表示为τ_{a λ} =βλ^－α，其中，α为气溶胶波长指数，β为混浊度系数.假设气溶胶程辐射与气溶胶光学厚度呈线性关系，则ρ_pa(λ)可以表达为(Aranuvachapun，1986)：

将阴影植被暗像元红、蓝波段ρ_pa(λ)代入式(8)，λ为中心波长，可快速求得Γ和δ.δ与波长指数α相当，对气溶胶粒度有指示作用，蕴含了气溶胶类型信息.由δ计算获得分子散射比例PER：

气溶胶单次散射反照率ω₀与分子散射比PER相关：

气溶胶散射相函数具有以下关系：

对于阴影植被像元，忽略太阳直射光，且由于τ_a较小，根据泰勒展开式将ρ_pa(λ)、ρ_dr(λ)、ρ_da(λ)中含τ_a的指数项按一次展开，带入式(2)，化简可得τ_a近似表达式：

则精确式可表达为：

实际应用时，先用式(12)算出光学厚度近似值τ_a^*，然后带入式(13)算出较精确值；再将得到的结果作为τ_a^*代入式(13)进行更精确计算，如此多次迭代计算可消除简化造成的误差，一般迭代两次后τ_a即可收敛.在缺少阴影植被的地区也可以利用浓密植被作为暗像元，同理可得浓密植被像元τ_a^*近似表达式：

精确式可表达为：

为了剥离气溶胶信息，气溶胶反演算法会对地表反射率进行近似表达.MODIS气溶胶产品算法利用浓密植被短波红外地表反射率对红、蓝波段地表反射率进行经验表达(Kaufman et al., 1988； Remer et al., 2005；Levy et al., 2007).ZY-3多光谱数据无法应用此类算法，因此，本文将植被像元的地表反射率表达为植被和裸土两个端元的线性光谱组合，通过NDVI解算植被占像元的比例，最终求解地表反射率(何颖清等，2010).此方法与BAER算法相似，但BAER算法直接将NDVI值作为植被端元丰度，计算可见光波段的地表反射率(von Hoyningen-Huene et al., 2003；2006).

暗像元作为区域气溶胶状况的模拟控制点，其选择的准确性将直接影响该区域气溶胶反演精度.由于大气状况分布不均，准确的暗像元越多、分布越广，对于大气状况空间分布的模拟就越细致越精确.本文采用先初选、后分区分规则优选的方法挑选阴影植被及浓密植被两类暗像元，以保证暗像元均匀密布于图像范围内.

阴影植被像元的NDVI值较同覆盖度的植被像元稍小，且红波段及近红外波段反射率均较小.则阴影植被像元的初判规则为：NDVI>0.3 且0.1<ρ_TOA(NIR)<0.25且ρ_TOA(RED)<0.10.近红外波段行星反射率大于0.1排除清洁和浑浊水体及混合像元，小于0.25排除浓密植被；红波段行星反射率小于0.10判断反射率低的暗像元(Richter et al., 2006).考虑到在气溶胶光学厚度达到2时，会使 NDVI 从 0.7 降低到 0.3 左右(王中挺等，2009)，因此，NDVI阈值设定为0.3.为了降低阈值设定宽松而错选混合像元的可能，增加周围8个邻域像元也必须同时满足该规则的判断，以此获得阴影植被像元初选集.一幅图像范围内大气状况不同，对NDVI的影响也不同，判断条件应当随区域而变.将图像分成N行×N列分区统计NDVI，去除小于NDVI均值的暗像元，此时获得的暗像元仍然是冗余而密集的，若将其全部带入插值计算，不单会造成算法效率低下，还会引入误差，需进一步细分区域筛选.假定M行×M列区域内大气状况不变，选取各小区域内最可靠的暗像元为代表.最优像元判断规则为红波段行星反射率最小，这是由于红波段受大气影响相对较小，各地物反射率在红波段越低，就越接近暗像元的假设.N取值应较大，才具有区域统计意义.而M取值应相对较小，才能确保暗像元具有区域代表性.默认情况下，N取值为500，M为20.浓密植被像元NDVI值较高，初判规则改为NDVI>0.4，其余选取过程与阴影植被选取一致.

为确保暗像元的准确性，最后对所有暗像元进行检验.若暗像元红、蓝波段的气溶胶波长指数α<0或α>4.08 或混浊度系数β<0，则剔除该像元.利用这些离散暗像元的τ_a值进行克里金插值就能获得整个图幅的气溶胶分布状况.

实验采用PCI Geomatica 9.0软件EASI模块编程实现气溶胶反演算法.对于香港地区采用海洋型气溶胶不对称因子参数(Gilabert et al., 1994).实验共选取了61636个暗像元，这些散布于不同区域、不同海拔的暗像元反映了子区内的大气状况，保证了大气参数的精确模拟.由ZY-3多光谱数据反演获得的550 nm处AOD呈现分布图(图 1a)可见，香港新界、九龙、湾仔繁华都市中心区AOD高于郊野山地区的趋势，且深圳市区的AOD明显高于香港地区.对比图 1a及图 1b可以发现，ZY-3反演的AOD结果与MODIS气溶胶产品趋势一致，均显示出深圳地区AOD显著高于香港地区的趋势.但分辨率为10 km的MODIS气溶胶产品显得较为粗糙，难以反映香港地区内部AOD的细节差异，甚至在某些区域由于缺少可靠的暗像元，出现结果缺失现象.由于ZY-3数据高分辨率的特点，成像时就能记录更多局地小型大气污染源信息，且在高地表反射率的城区包含更多可靠的浓密植被“纯净”暗像元，因此，ZY-3数据对城市区AOD的反演更精细且完整.以蓝框(图 1c)及红框(图 1d)子区为例，新界元朗地区、九龙及香港岛北部湾仔区均为路网、楼宇密集区，由于交通及生活污染排放，其AOD明显高于周边植被覆盖良好的郊野山地公园.说明ZY-3反演结果在揭示繁华都市区内部的AOD差异，识别城市内部污染源方面更具优势.另外值得一提的是，AOD值较高的区域大都集中于100 m高程以下，其原因一方面是由于地势低的城区较山地区大气污染排放量大，另一方面是由于山地对污染物的扩散具有阻挡作用，且气溶胶自身具有随高度增大呈指数衰减的规律.因此，污染物易在近地面聚集，这与Bilal(2013)及黄文声等(2009)的研究结论相似.

图 1(Fig. 1)

图 1 ZY-3气溶胶光学厚度与MODIS气溶胶产品MOD04(C051)对比图(a.ZY-3 数据550 nm处AOD;b.MODIS气溶胶产品MOD04(C051);c.蓝色子区AOD;d.红色子区AOD) Fig. 1 Comparisons of AOD image of ZY-3 MUX and MODIS aerosol product(MOD04 C051)(a.AOD image of ZY-3 MUX，b. MODIS aerosol product(MOD04 C051)，c. subset AOD of blue frame，d. subset AOD of red frame)

计算获得气溶胶光学参数后即可带入式(2)反算地表真实反射率信息，去除大气影响.受到大气散射的影响，纠正前图像蓝光波段反射率极高(图 2a)，纠正后图像蓝波段的反射率明显降低，蓝灰色消除(图 2b).纠正后图像恢复地物原貌，对比度及清晰度增强.取红框内新界元朗区域进行对比，都市区及山地区大气纠正的效果均良好，纠正后建筑物、道路轮廓清晰，河涌及池塘内水体细节信息得到恢复.

图 2(Fig. 2)

图 2 大气纠正前后对比图(R:Red，G:Green，B: blue)(a.大气纠正前，b.大气纠正后，c.子区大气纠正前，d.子区大气纠正后) Fig. 2 Comparisons of the true-color composites of the original and corrected images(a.overall original，b.overall corrected，c.subset original，d.subset corrected)

统计图像大气纠正前后较清洁水体和浓密植被像元各50个，考察均值光谱特征，可间接评价方法的有效性，结果如图 3所示.纠正后，水体和植被在可见光波段的反射率都降低了，恢复低反射率特征.由于大气散射随波段递减，因此，减幅最大为蓝波段.去除大气影响后，植被光谱显现出绿光反射峰特性.由于水体在红外波段的强吸收性，纠正后反射率急剧降低，而植被反射率在红外波段反而略有提升.其原因主要是大气散射增大了低反射率地物如水体的行星反射率，而对于高反射率的植被受大气散射、吸收的削弱作用更明显，降低了行星反射率，因此，纠正后反射率反而增大了.

图 3(Fig. 3)

图 3 大气纠正前后反射率比较(a.浓密植被，b.清洁水体; 波段1：蓝波段，波段2：绿波段，波段3：红波段，波段4：近红外) Fig. 3 Reflectance derived from with and without atmospheric correction(a. dense vegetation，b. clean water，B and 1：Blue，B and 2：Green，B and 3：Red，B and 4：NIR)

1)本研究方法利用红、蓝波段阴坡植被暗像元程辐射信息，估算气溶胶相函数和散射比，能更客观地从图像中获取气溶胶光学参数，并考虑了参数的空间异质性.计算暗像元气溶胶光学厚度时，本方法考虑了地表对天空漫射光的反射贡献，通过NDVI按线性端元分解的方法估算地表反射率，相较Gilabert模型更具完备性.

2)香港地区ZY-3数据AOD反演结果与MODIS气溶胶C051产品总体趋势一致，但ZY-3数据能够反映繁华都市区内部的AOD差异，识别城市内部大气污染源.

3)基于图像自身反演的AOD进行ZY-3数据大气纠正，从纠正后图像视觉效果增强、统计水体及浓密植被光谱特征与先验知识相符这两方面证明了算法的合理性.

4)算法对辐射传输过程进行简化，仅基于图像就能自动实现气溶胶光学厚度反演及大气纠正.然而算法只考虑了一次散射，忽略了周围像元的临近效应，在地形效应及气溶胶模式等方面仍具有改进空间.另外，此次实验仅利用了一个时相的ZY-3数据反演结果与MODIS气溶胶C051产品进行对比，今后应进一步开展完备的验证工作.