一种协同反演气溶胶与水汽含量的高光谱图像大气校正算法

佃袁勇，方圣辉，徐永荣

(1．华中农业大学园艺林学学院，武汉 430070;2．武汉大学遥感信息工程学院，武汉 430079)

0 引言

大气中的气溶胶和水汽是影响航空与航天遥感图像数据质量的重要因素之一［1－5］。大气校正的目的就是要消除大气中各种成分对遥感信号的影响，还原真实的地表辐射信息。大气的状况瞬息万变，很难同步测量卫星过境时大气的参数;而如果能直接从遥感信息中反演出大气信息，就能更好地进行大气校正。大气中的 O3，O2，CO2，NO2和 CH4等含量相对稳定，只有气溶胶和水汽含量变动较大，因而大气校正的关键就在于估算气溶胶模式、气溶胶光学厚度(aerosol optical thickness，AOT)与水汽含量(water vapor content，WV)［3－9］。目前大气校正方法主要有基于图像特征的相对订正法、基于地面线性回归模型法、大气辐射传输模型法和复合模型法等［7－17］。其中，大气辐射传输模型法模拟太阳辐射信号通过大气并与地表相互作用后回到传感器的信号，能较精确地描述水汽和气溶胶的作用［8－12］，因此被广泛应用于遥感图像的大气校正。但现阶段基于大气辐射传输模型的大气校正算法对AOT和WV是分开进行反演的——先利用水汽吸收波段估算 WV［18－20］，在消除水汽影响后再考虑 AOT 的反演;而且在考虑气溶胶的影响时，假设气溶胶为单一的模式，没有考虑气溶胶模式的差异问题。但实际中气溶胶模式是变化的，对气溶胶模式的选择是最大的误差源［3－5，10，21］。现有研究成果表明，大陆型和城镇型气溶胶模式在蓝波段和红波段的反射率最大差别可达到 13%［22－24］。

本文在前人研究基础上，提出了一种协同反演AOT与WV的高光谱遥感图像大气校正算法。该算法以6S辐射传输模型［25］(Version 4．1)构建查找表，在估算气溶胶模式的基础上，采用循环迭代方式协同反演AOT与WV;并以EO－1卫星Hyperion高光谱图像为例，验证该算法的有效性。

1 实验数据

以武汉市Hyperion高光谱数据为数据源，验证本文提出的大气校正算法。Hyperion是地球观测卫星EO－1上搭载的高光谱传感器，共有242个波段(其中44个波段没有定标，有效波段只有198个);其中第8—57波段采用可见光的定标系数，第77—224波段采用短波红外的定标系数;空间分辨率为30 m，光谱分辨率为10 nm。本文使用的数据是2004年4月2日上午10:40经过武汉市上空获取的Hyperion卫星数据(行列号为123/39)。

Hyperion数据信号转化为表观反射率的公式为

式中:ρλ为λ波段的表观反射率;L为光谱辐亮度;d为日地归一化距离;E为大气上界太阳光辐照度;θ为太阳天顶角;DN为图像像元值;S为定标系数，对于可见光波段，S=40;对于短波红外波段，S=80。

鉴于没有地面同步实测数据可供对比，本文选择与EO－1卫星Hyperion准同步过境的TERRA卫星MODIS的AOT与WV产品，用于评估本文提出的反演算法。由于MODIS的AOT产品的分辨率为10 km，WV产品的分辨率为1 km(选择由近红外波段反演的WV产品)，其分辨率远小于Hyperion高光谱图像的分辨率;因此，本文根据MODIS产品提供的经纬度坐标，选取在Hyperion图像范围内的AOT和WV的均值作为验证数据(表1)。

表1 Hyperion图像范围内从MODIS产品中提取的AOT与WV数据Tab．1 AOT and WV data extracted from MODIS products in Hyperion imange extent

2 研究方法

根据6S辐射传输理论，卫星接收到的辐射信号包含大气直接散射、地表目标反射和地表目标周围背景反射部分。假设地表为朗伯面，则卫星入瞳处信号(即表观反射率)ρTOA为

式中:ρTOA为传感器接收到的单个像元的反射率;Tg(O3，O2，CO2，NO2，CH4)为 气 体 分 子 引 起 的 大 气 透 过 率;ρR+A为大气瑞利和分子散射;TR+A为瑞利散射和气溶胶引起的大气透过率;Tg(H2O)为水汽引起的大气透过率;S为大气半球反射率;ρs为像元地表反射率。所有的变量均与波长相关。

由ρTOA反演地表反射率ρs的过程需要首先估算大气参数，即 Tg(O3，O2，CO2，NO2，CH4)，ρR+A，S，TR+A和T g(H2O)。将 Tg(O3、O2，CO2，NO2，CH4)和 ρR+A归纳为一个变量 ρpath，将 T g(O3、O2，CO2，NO2，CH4)，T g(H2O)和 T R+A归纳为一个变量T，则式(2)可简化为

式中:ρpath，T和 S分别为与大气廓线，水汽含量(WV)和气溶胶模式、气溶胶光学厚度(τ)有关的变量。当给定大气廓线、气溶胶模式、WV和τ时，利用6S模型即可计算ρpath，T和S这3个变量;进而在给定星上反射率后，即可得到地表反射率。于是大气校正的问题转化为:在已知星上反射率的情况下，求解大气廓线、气溶胶模式、WV和AOT参数，亦即求解ρpath，T和S。其中，大气廓线可根据遥感数据获取时间以及所在地区，选择6S模型中已有的廓线数据。因此，WV、气溶胶模式和AOT是大气校正需要反演的3个参数。本文提出了一种基于AOT和WV协同反演的大气校正算法，该算法的流程如图1所示。

图1 大气校正算法流程图Fig．1 Flow chart of atmospheric correction algorithm

算法的关键步骤包括气溶胶模式反演、查找表生成、WV反演以及暗目标和蓝波段(0．48μm)及红波段(0．66μm)地表反射率确定等。

2．1 反演气溶胶模式

以6S模型中定义的4种类型的气溶胶(即沙尘类型(dust－like，DL)、海洋型(oceanic，OC)、水溶型(water－soluble，WS)和煤烟类(soot，SO))为基础，将这4种类型气溶胶按一定比例组合，建立不同的气溶胶模式;结合MODIS提供的AOT与星上反射率数据，反演气溶胶模式。根据武汉市所处的地理位置和环境等情况，武汉市城区气溶胶的组成在城市型气溶胶(体积比为沙尘型占17%，水溶型占61%，煤烟类占22%)和大陆型气溶胶(体积比为沙尘型占70%，水溶型占29%，煤烟类占1%)之间，因3种组分之和应等于100%，故3个变量并不是互相独立的。首先把沙尘型和烟尘型组分各分为10个等级，水溶型组分可根据前两者含量计算得到;然后根据MODIS提供的AOT，利用这100种气溶胶模式，根据查找表分别反演地面暗目标在红波段和蓝波段的星上反射率和。假定550 nm处的光学厚度为τ，卫星观测的表观反射率在红波段和蓝波段分别为ρR和ρB，当满足条件时(式中ε为表观反射率的误差允许范围)，得到的气溶胶模式就是当前大气的气溶胶模式。

2．2 反演水汽含量

在近红外通道，当忽略地表反照率的影响，将其他的吸收特性与瑞利、气溶胶特性融合在一起考虑时，则辐射传输方程(2)可简化为

式中:T为大气中水汽含量导致的透过率T(H2O)与其他成分的综合透过率Tz的乘积;ρs为地表反射率值。

而水汽通道的透过率可以表示为与水汽含量的关系［18］，其模型为

式中:ω为水汽含量;α和β为常量参数。

利用6S和MODTRAN等模型可以很好地建立透过率与水汽含量的关系，因此求解的关键是如何求透过率。现有的水汽反演算法中，将式(5)中的ρpath项忽略掉，同时认为在0．85～1．25 μm 之间的地面反射率满足线性关系［18］，且水汽以外窗口的大气透过率为1;将大气透过率表达为2个波段或3个波段的比值，进而利用水汽查找表求得水汽含量［19］。而在实际过程中，由于受气溶胶的厚度和类型等影响，ρpath并不为0，且水汽以外窗口的大气透过率也不是1，因此利用该算法反演的水汽含量存在一定误差。

本文在以上研究基础上，考虑了气溶胶的影响因素，对3个波段的水汽反演算法进行了改进。针对Hyperion高光谱的卫星数据中940 nm和1 140 nm处水汽吸收带具有的多个波段，根据赵祥等［7］的建议，选择Hyperion第80波段(942．73 nm)为水汽吸收波段，选择第52波段(874．53 nm)和第110波段(1 245．36 nm)为水汽弱吸收波段;采用上述3个波段比值模型，具体构建水汽透过率T(H2O)的计算公式，即

式中:ρTOA(80)，ρTOA(52)和 ρTOA(110)分别为 Hyperion 第80，52 和110 波段的星上反射率;ρpath(80)，ρpath(52)和ρpath(110)分别为Hyperion第80，52和110波段的大气程辐射;C1=0．82;C2=0．18。

2．3 生成查找表

在确定研究区域的大气气溶胶模式后，利用6S模型，根据相应的观测几何参数，可以建立WV，τ与ρpath，T，S之间对应关系的查找表。在构建查表时，WV 的变化范围在 0．2 ～4．2 g/cm2，以0．2 g/cm2递增;τ的变化区间在0．05～5．0，其中 τ＜2时以0．2递增，τ＞2时以0．5递增。针对Hyperion图像，在给定的观测几何参数下，对每一个τ和WV值，假定3 个地表反射率值(ρs=0．0，0．5，0．8)，用 6S 模型进行运算，可得到3个星上反射率。将3组对应的地表反射率和星上反射率代入式(3)，即可求得ρpath，T和S。查找表即可表示为WV和τ与地表反射率和星上反射率的关系表。这样，对于大气参数(AOT和WV)、地表反射率和星上反射率这3组变量，只要知道了其中2组，便可以推导出第3组变量。

2．4 确定暗目标

在高光谱图像大气校正中，图像中的暗目标能用于探测AOT［5］。对于绿色植被、黑色土壤等暗目标区域在红波段(0．66μm)和蓝波段(0．48μm)的地表反射率，根据Levy等［5］的建议，在暗目标地区，红波段)与蓝波段)地表反射率与2．12 μm处的暗目标反射率)之间的关系可以表示为NDVISWIR和散射角Θ的函数，其函数关系为

其中，

式(10)—(11)中:

式(15)—(16)中:ρ1．24和 ρ2．12分别为 1．24 μm 和2．12μm处的星上反射率;θ0为太阳天顶角;θ为传感器观测天顶角;φ为太阳与传感器间相对方位角。

本文在确定暗目标及红波段和蓝波段的地表反射率时，直接采用了上述方法。

2．5 过程循环迭代反演

在确定了研究区域气溶胶模式并利用6S模型建立了查找表后，采取2．4节中的方法，利用高光谱图像数据确定暗目标，估算蓝波段和红波段的地表反射率;然后采用循环迭代的方式反演AOT与WV。具体过程如下:

1)给定初始光学厚度τ0，利用2．2节中的方法计算水汽含量初值WV0。

2)根据暗目标的星上反射率及估算的地表反射率(见2．4节)，利用已经建立的查找表(见2．3节)，得到气溶胶光学厚度τi。

3)根据τi，结合2．2节中水汽含量反演方法，重新计算水汽含量值WVi。

4)更新当前气溶胶光学厚度τi+1及水汽含量值WVi+1。

5)比较2次获得的AOT及WV值，如果2次的差值在误差允许范围内，则进入步骤6);否则，返回步骤2)。

6)输出WV与AOT，并计算每个像元的地表反射率。

3 结果与分析

3．1 气溶胶模式反演结果

以同一天获取的MODIS原始数据及气溶胶产品，获得研究区域的红波段和绿波段地表反射率、星上反射率以及相对应的AOT。利用气溶胶模式反演算法，反演得到武汉市的气溶胶模式为:沙尘型所占比例为17%，水溶型所占比例为82%，煤烟类所占比例为1%;此时式(4)中的ε为0．004。该气溶胶模式介于大陆型和城市型气溶胶之间，与武汉市的地理条件和环境(湿度较大，沙尘、煤烟较少)比较相符。

3．2 AOT和WV的验证与比较

利用暗目标确定方法，在整个研究区范围内共确定2 196个暗目标像元;根据这些暗目标像元的地表反射率与星上反射率数据，采用本文提出的算法反演大气参数。为比较不同算法的反演效果，利用ENVI软件中FLAASH大气校正模块对同一Hyperion高光谱图像进行大气校正，以比较用FLAASH反演的AOT和WV与本文算法反演的结果。同时，以同一天获取的MODIS的AOT和WV产品值为标准值(因没有当天的地面实测数据，故无法与地面数据对比)计算相对误差，其结果如表2所示。

表2 AOT与WV反演结果Tab．2 Inversion results of AOT and WV

从表2可以看出，与MODIS的值相比本文方法对WV的反演精度，相对误差在3．0% ～37．6%之间，平均相对误差20．3%，精度略高于FLAASH算法的反演结果;对AOT的反演精度，相对误差在5．3% ～35．5%，平均相对误差 11．8%，精度略低于FLAASH算法的反演结果。

3．3 大气校正前后光谱反射率比较

图2显示了大气校正前、后的效果，其中图2(a)为大气校正前的Hyperion高光谱图像，图2(b)与(c)分别为FLAASH算法与本文算法的大气校正结果。

图2 大气校正前后效果对比Fig．2 Comparsion between results before and after atmospheric correction

图3显示了大气校正前、后Hyperion高光谱图像 中典型植被、水体和水泥路面光谱反射率差异。

图3 大气校正前后典型地物的光谱反射率Fig．3 Spectrum reflectance of typical objects before and after of atmospheric correction

对比图3中大气校正前、后的效果可以看出，气溶胶对可见光波段(400～700 nm)，特别是对蓝波段影响较大;水汽的影响主要在近红外波段，特别是在820 nm，940 nm和1 135 nm附近。经过大气校正后，可基本消除气溶胶及水汽的影响。对比本文算法与FLAASH算法大气校正后典型地物的光谱反射率曲线可以看出，植被的反射率相差最大为0．09，水体的反射率相差最大为0．07，水泥路面反射率相差最大为0．05，表明本文算法的大气校正效果与FLAASH算法的相当。

3．4 循环迭代效果比较

本文采用循环迭代的思路协同反演AOT与WV，对整个研究区域内的2 196个暗目标像元计算了相邻2次迭代的误差。图4显示了迭代次数与相邻2次迭代的AOT和WV的平均误差关系。

图4 反演AOT和WV的迭代次数与相邻2次迭代误差的关系Fig．4 Relationship between iteration times and adjacent two iterative error in retrieval AOT and WV

从图4可知，随着迭代次数的增加，相邻2次迭代的误差逐渐减小，一般迭代3～4次即可收敛，此时AOT和WV的误差均＜0．05。但是，每个像元都要循环迭代3～4次，会增大数据处理的运算量。

4 结论

1)大气校正是高光谱图像定量反演地表参数的前提。本文提出了一种协同反演大气气溶胶光学厚度(AOT)与水汽含量(WV)的方法，充分利用高光谱数据本身的光谱特点，在同时考虑气溶胶模式、AOT和WV这3个因素的综合影响基础上，采用循环迭代的思想，基于6S辐射传输模型，反演大气参数及地表反射率，弥补了现有反演算法中没有同时考虑AOT与WV的不足。

2)从对武汉市Hyperion高光谱图像大气校正的效果来看，本文提出的算法能较好地校正大气中气溶胶与水汽对高光谱图像的影响，且反演过程中所有的输入均来自图像数据本身或6S辐射传输模型，无需输入额外的参数。

由于本文提出的算法采用了循环迭代的思路，运算量较大，增加了大气校正的时间。同时，因缺少地面的实地验证数据，仅将本文方法的大气校正结果与FLAASH大气校正结果和MODIS数据产品进行了对比。在今后的工作中，需要进一步优化算法，提高算法的效率，并进行同步的地面验证工作。

［1］ Barducci A，Guzzi D，Marcoionni P，et al．Algorithm for the retrieval of columnar water vapor from hyper spectral remotely sensed data［J］．Applied Optics，2004，43(29):5552－5563．

［2］ Li BG，Yuan H S，Feng N，etal．Comparing MODIS and AERONET aerosol optical depth over China［J］．International Journal of Remote Sensing，2009，30(24):6519－6529．

［3］ Levy R C，Remer L A，Dubovik O．Global aerosol optical properties and application to moderate resolution imaging spectroradiometer aerosol retrieval over land［J］．Journal of Geophysical Research:Atmospheres，2007，112(D13):D13210．

［4］ 夏 双，阮仁宗，张 月，等．气溶胶光学厚度对蓝藻水华信息提取的影响［J］．国土资源遥感，2013，25(1):33－38．doi:10．6046/gtzyyg．2013．01．06．Xia S，Ruan R Z，Zhang Y，et al．Effects of aerosol optical thickness on extracting cyanbacteria bloom［J］．Remote Sensing for Land and Resources，2013，25(1):33－38．doi:10．6046/gtzyyg．2013．01．06．

［5］ Levy R C，Remer L A，Mattoo S，et al．Second－generation operational algorithm:Retrieval of aerosol properties over land from inversion of moderate resolution imaging spectroradiometer spectral reflectance［J］．Journal of Geophysical Research:Atmospheres，2007，112(D13):D13211．

［6］ Remer L A，Kaufman Y J，TanréD，et al．The MODIS aerosol algorithm，products，and validation［J］．Journal of the Atmospheric Sciences，2005，62(4):947－973．

［7］ 徐永明，覃志豪，陈爱军．基于查找表的MODIS逐像元大气校正方法研究［J］．武汉大学学报:信息科学版，2010，35(8):959－962．Xu Y M，Qin ZH，Chen A J．A pixel－by－pixel atmospheric correction algorithm for MODIS data based on look－up table［J］．Geomatics and Information Science of Wuhan University，2010，35(8):959－962．

［8］ 赵 祥，梁顺林，刘素红，等．高光谱遥感数据的改正暗目标大气校正方法研究［J］．中国科学:地球科学，2007，37(12):1653－1659．Zhao X，Liang S L，Liu SH，et al．Atmospheric correction with Hyperspectral image data based on modified dark object method［J］．Scientia Sinica Terrae，2007，37(12):1653－1659．

［9］ 郑求根，权文婷．基于暗像元的Hyperion高光谱影像大气校正［J］．光谱学与光谱分析，2010，30(10):2710－2713．Zheng Q G，Quan W T．Application of dark pixels atmospheric correction algorithm to Hyperion imageries［J］．Transactions of Spectroscopy and Spectral Analysis，2010，30(10):2710－2713．

［10］ 胡方超，王振会，张 兵，等．遥感试验数据确定大气气溶胶类型的方法研究［J］．中国激光，2009，36(2):312－317．Hu FC，Wang ZH，Zhang B，et al．Study on method for determining atmospheric aerosol type using remote sensing experimental data［J］．Chinese Journal of Lasers，2009，36(2):312－317．

［11］ 唐洪钊，晏 磊，李成才，等．基于MODIS高分辨率气溶胶反演的ETM+影像大气校正［J］．地理与地理信息科学，2010，26(4):12－15．Tang H Z，Yan L，Li C C，et al．Atmospheric correction for ETM+imagery based on high resolution aerosol optical depth retrieved from MODIS data［J］．Transactions of Geography and Geo－Information Science，2010，26(4):12－15．

［12］ 袁金国，牛 铮，王锡平．基于FLAASH的Hyperion高光谱影像大气校正［J］．光谱学与光谱分析，2009，29(5):1181－1185．Yuan JG，Niu Z，Wang X P．Atmospheric correction of Hyperion hyper spectral image based on FLAASH［J］．Transactions of Spectroscopy and Spectral Analysis，2009，29(5):1181－1185．

［13］ Liang SL，Fang H L．An improved atmospheric correction algorithm for hyper spectral remotely sensed imagery［J］．Geoscience and Remote Sensing Letters，IEEE，2004，1(2):112－117．

［14］ 程 晨，陈 健，李鑫慧．基于TM图像的南京市气溶胶光学厚度反演［J］．国土资源遥感，2013，25(3):90－96．doi:10．6046/gtzyyg.2013．03．16．Cheng C，Chen J，Li X H．Retrieving aerosol optical depth over Nanjing City based on TM image［J］．Remote Sensing for Land and Resources，2013，25(3):90－96．doi:10．6046/gtzyyg．2013．03．16．

［15］ 杨 磊，潘志强，傅俏燕，等．资源三号卫星多光谱数据的大气校正研究［J］．光谱学与光谱分析，2013，33(7):1903－1907．Yang L，Pan ZQ，Fu Q Y，etal．Research on the atmospheric correction for ZY－3MUX image［J］．Spectroscopy and Spectral Analysis，2013，33(7):1903－1907．

［16］ 方 莉，余 涛，顾行发，等．北京地区HJ－1卫星CCD数据的气溶胶反演及在大气校正中的应用［J］．遥感学报，2013，17(1):151－164．Fang L，Yu T，Gu X F，et al．Aerosol retrieval and atmospheric correction of HJ－1 CCD data over Beijing［J］．Journal of Remote Sensing，2013，17(1):151－164．

［17］ 杨贵军，黄文江，刘三超，等．环境减灾卫星高光谱数据大气校正模型及验证［J］．北京大学学报:自然科学版，2010，46(5):821－828．Yang G J，HuangW J，Liu SC，etal．Research on modeling and validating of atmospheric correction for HJ－1A hyper spectral imager data［J］．Transactions of Peking University:Natural Science Edition，2010，46(5):821－828．

［18］ Kaufman Y J，Gao BC．Remote sensing of water vapor in the near IR from EOS/MODIS［J］．IEEE Transactions on Geoscience and Remote Sensing，1992，30(5):871－884．

［19］ Gao B C，Kaufman Y J．Water vapor retrievals using moderate resolution imaging spectroradiometer(MODIS)near－infrared channels［J］．Journal of Geophysical Research:Atmospheres，2003，108(D13):4389．

［20］ 姜立鹏，覃志豪，谢 雯．针对MODIS近红外数据反演大气水汽含量研究［J］．国土资源遥感，2006，18(3):5－9．doi:10．6046/gtzyyg．2006．03．02．Jiang L P，Qin ZH，XieW．Retrieving atmospheric water vapor from MODIS near infrared data［J］．Remote Sensing for Land and Resources，2006，18(3):5－9．doi:10．6046/gtzyyg．2006．03．02．

［21］ 周春艳，柳钦火，唐 勇．MODIS气溶胶C004、C005产品的对比分析及其在中国北方地区的适用性评价［J］．遥感学报，2009，13(5):863－872．Zhou C Y，Liu Q H，Tang Y．Comparison between MODIS aerosol product C004 and C005 and evaluation of their applicability in the north of China［J］．Transactions of Journal of Remote Sensing，2009，13(5):863－872．

［22］ 佃袁勇．高光谱数据反演植被信息的研究［D］．武汉:武汉大学，2011．Dian Y Y．Inversion Vegetation Information Based on Hyper spectral Remote Sensing Data［D］．Wuhan:Wuhan University，2011．

［23］ Laurent V CE，VerhoefW，Clevers JG PW，etal．Estimating forest variables from top－of－atmosphere radiance satellite measurements using coupled radiative transfer models［J］．Remote Sensing of Environment，2011，115(4):1043－1052．

［24］ 陈新芳，陈镜明，安树青，等．不同大气校正方法对森林叶面积指数遥感估算影响的比较［J］．生态学杂志，2006，25(7):769－773．Chen X F，Chen JM，An SQ，et al．Comparison of different atmospheric correction models in their effects on Landsat TM estimation of forest leaf area index［J］．Chinese Journal of Ecology，2006，25(7):769－773．

［25］ Vermote E F，TanréD，DeuzéJL，et al．Second simulation of the satellite signal in the solar spectrum，6S:An overview［J］．IEEE Transactions on Geoscience and Remote Sensing，1997，35(3):675－686．