| 标题 | 基于MODTRAN的高光谱快速大气校正方法 |
| 范文 | 郭飞 摘 要： 高光谱遥感成像过程中一般会受大气成分的影响，为了利用高光谱数据进行有效的地物识别，这些影响必须消除。详细介绍大气校正的一个过程：首先采用大气预处理微分吸收法（APDA）和暗像元法，通过并行调用MODTRAN先后建立水汽和能见度的查找表，从而计算出图像中水汽含量和能见度，并利用得出的这些参数建立反射率查找表，最终对高光谱图像进行逐像元校正。该过程结合辐射传输模型法和基于图像信息的校正方法的优点，不需要实时的大气参数也能进行绝对大气校正。同时，加入并行处理使运行速度显著提高。通过对校正后的数据进行对比分析，结果表明校正后的光谱曲线与地面实物和其他软件校正后的光谱特征相吻合，该方法能够满足大气辐射校正的基本要求。 关键词： 大气校正； MODTRAN； APDA； 暗像元法 中图分类号： TN911?34 文献标识码： A 文章编号： 1004?373X（2015）10?0023?05 0 引 言 20世纪80年代初期成像光谱（高光谱）技术的出现，使光学遥感进入了崭新的阶段。高光谱遥感能在电磁波谱的可见光和红外区域， 获取许多非常窄的光谱图像数据，为每个像元提供数十至数百个窄波段（一般小于10 nm）的光谱信息，产生一条完整而连续的光谱曲线，以便于识别出地物特征，具有波段多、光谱分辨率高、图谱合一等特点。大气校正是遥感信息定量化过程中不可缺少的环节。在高光谱图像的拍摄过程中，由于大气的存在，进入大气的太阳辐射会发生反射、折射、吸收、散射[1]，其中对传感器接受影响较大的是吸收和散射。由于大气的吸收和散射，使得进入到传感器中的信号减弱[2]。同时，大气的散射光也有一部分直接或经过背景地物反射进入到传感器时却又使得原型号增强了。为了消除这些因素的影响，必须对高光谱数据进行大气辐射校正。 目前大气校正的方法主要有基于图像特征的经验线性法、相对校正方法和基于物理的辐射传输模型法[3]。经验线性法需要做大量的地物反射率测量工作，通过统计的方法建立传感器接收到的辐亮度与地物反射率的线性回归方程。基于图像特征的相对校正方法得到的不是地物的真实反射率，而是相对反射率。它不需要实际地面光谱和大气参数的测量，可以应用于历史数据。它校正后的结果不是具体的反射率或辐亮度，相同的结果只表示相同的反射率，至于多少不能给定[4]。辐射传输模型法建立在大气辐射传输方程的基础上，在进行大气辐射校正时需要输入各种大气参数（例如：卫星过境时的地物反射率、大气的能见度、太阳天顶角和卫星传感器的标定参数等），获得这些参数的代价昂贵，又不易实现，尤其对于历史存档数据，获取这些参数更是不可能的[5]。陶东兴，赵慧洁通过调用MODTRAN分别建立水汽吸收通道和2个参考通道的2个查找表[6]，搜索这2个查找表使用三次样条插值方法计算水汽含量，实验结果表明生成的水汽含量分布图精度与MODTRAN计算结果只有0.1%的误差。另外还有暗像元法，暗目标方法通过确定暗目标在可见光波段的地表反射率来估计气溶胶参数。赵祥等选取小麦作为新的暗目标[4]，使用波段间的线性回归与插值的方法对常规暗目标方法进行改进，从而将暗目标方法扩展应用于高光谱遥感数据的大气校正。彭妮娜等提出了浓密植物像元660 nm与830 nm波段处的线性关系进行气溶胶含量反演[7]，为获取气溶胶含量提供了一种新途径。 1 概 述 为了保证大气校正的精确性和程序运行速度，结合前人的工作经验，总结出一套大气辐射校正方法，形成一整套大气辐射校正的具体流程，在流程中加入并行计算提高运行效率。 该方法具体思路是通过大气预处理微分吸收法（APDA法）并调用MODTRAN建立的水汽查找表来获得水汽含量，然后采用暗像元法并调用MODTRAN建立的能见度查找表获得能见度，最后利用水汽含量和能见度参数调用MODTRAN建立反射率查找表，计算FLAASH模型中的大气辐射校正的参数从而实现大气辐射校正的目标。 运用该方法开发了用于大气辐射校正的软件。并对调用MODTRAN建立查找表进行并行处理，提高运行效率。该软件包括水汽反演、能见度反演、大气辐射校正这三大模块。为了验证这三大模块的精确性，对东天山地区高光谱数据（已完成辐射定标和几何校正）进行大气辐射校正，将校正后的光谱和地面采集的光谱进行对应像元对比，分析校正结果的准确性。本文方法的大体思路流程见图1。 图1 大气校正流程图 2 MODTRAN快速查找表并行处理 2.1 调用MODTRAN计算大气参数 在调用MODTRAN后会生成tape7文件，该文件得到一系列参数（FREQ，TOTAL RAD，GRND RFLT，DRCT RFLT等）。 利用这些参数可以计算如下的大气校正所需的参数： [λ=10×107FREQ] （1） [Lp=（TOTAL RAD）0] （2） [S=（（GRND RFLT）1.0-（GRND RFLT）0）（（GRND RFLT）1.0-（GRND RFLT）0）-→ ←-2?（（GRND RFLT）1.0-（GRND RFLT）0）（（GRND RFLT）1.0-（GRND RFLT）0）] （3） [Gb=（（TOTAL RAD）1.0-（TOTAL RAD）0）·（1-S）] （4） [Gt=（（DRCT RFLT）1.0-（DRCT RFLT）0）·（1-S）] （5） 式中：[λ]表示波长；[Lp]表示程辐射；[S]为大气球面反照率；[Gb]，[Gt]为辐射校正参数。 2.2 并行处理 并行处理是一个对数据进行分片很有效的方法，在调用MODTRAN时可以将tape5文件进行分片，独立生成tape7文件。 以建立水汽查找表为例，为了提高运行速度，在写tape5文件时，将12个不同的水汽含量和3种不同的反射率共36种情形分别写入到12个文件夹中，在这12个文件夹中独立调用以达到减少程序运行时间的目的，在建立查找表需要整个数据时再将这些数据合并。而且高光谱数据处理的典型特点是像素间、波段间依赖关系比较小，相互独立性强，在不同波段数据间也可以进行分片。 因此在3个模块中，涉及了逐个波段计算时，一般都对每个波段进行并行计算处理。以建立水汽查找表为例，并行操作的流程图如图2所示。 3 水汽含量反演 水汽含量的影响包括对大气透过率的影响和程辐射的影响，但是它却不是惟一的因素，气溶胶和大气分子的散射也对大气透过率有不可忽视的影响。本系统为了区分这两种因素的影响，根据水汽吸收具有波段选择性，而气溶胶和大气分子的散射是关于波长的缓变函数这一特点，选用大气预处理微分吸收法（Atmospheric Precorrected Differential Absorption，APDA）[7]。 3.1 大气预处理微分吸收法 该方法利用1 135 nm处的水汽吸收波段及其邻近的非水汽吸收波段的比值来获取大气水汽柱含量。在特定的大气条件下，水汽含量值与参数[RAPDA]有如下经验关系： [RAPDA=e-（γ+α（PW）β）] （6） 式中：PW是水汽含量；[α，β，γ]是未知的参数，在特定不变的大气条件下[α，β，γ]的值不变；[RAPDA]的计算方法如下： [RAPDA=[Lm-Lp，m]LIR（[λr]j，[Lr-Lp，r]j）λm] （7） 式中：[Lm]是水汽吸收通道的辐射值；[Lp，m]是水汽吸收通道的程辐射，在MODTRAN中通过设置地面反射率为0来获取各波段的程辐射；[Lr]是水汽参考通道的辐射值；[Lp，r]是水汽参考通道的程辐射；[LIR（[λr]j，[Lr-Latm，r]j）λm]是指通过线性拟合参考通道的辐射值，估计在水汽吸收通道应有的辐射值。 图2 建立水汽查找表流程图 3.2 水汽初始含量选择 水汽含量的初始值的选择直接关系着最终结果的正确性。关于如何设置水汽含量初值，ENVI的帮助文档中给出了设置方法：如果没可用的水汽信息，则根据已知的或预期的大气温度来设置水汽含量。如果温度是未知的，则根据数据获取时间的季节性来设置。总的来说，它是基于一个季节性纬度表面温度模型的设置方法[8]。 3.3 水汽含量反演流程图 利用大气预处理微分吸收法求解水汽含量具体流程图如图3所示。 图3 水汽含量反演流程图 4 能见度反演 大气中的气溶胶粒子将会吸收和散射太阳辐射，从而使卫星传感器接收到的辐射信号产生衰减，是影响能见度的主要因素。本项目基于精确度需要以及避免地面辐射定标的巨大工作量，选用暗像元法反演能见度。所谓暗像元是指某些波段反射率接近为0的像元。暗像元可以选择水体、密度较高的植被森林或阴影地区。 4.1 暗像元的选取 该方法的关键点在于暗像元的选取，暗像元选择的科学、正确与否关系到算法的精度。Kaufman等首先采用归一化植被指数[9]（Normalized Difference Vegetation Index，NDVI）进行暗目标确认： [NDVI=ρ860-ρ660ρ860+ρ660>0.7ρ2100<0.05] （8） 式中：[ρ860，ρ660，ρ2100]分别表示在860 nm波段、660 nm波段、2 100 nm波段的地物反射率，若满足式（8）则是暗像元，否则不是暗像元。 但是由于NDVI对气溶胶的变化较为敏感[10]，采用对大气不敏感的大气耐抗植被指数[11]（Atmospherically Resistant Vegetation Index，ARVI），在已得到的暗像元位置上继续进行再次选择得出最终的暗像元位置。该指数是通过地面调查得到地物类别在影像中的精确位置，经过分析后确定暗目标选取方式的经验公式如下： [ARVI=ρ*884+ρ*467-2ρ*660ρ*884-ρ*467+2ρ*660>0.5ρ*2042<0.04ρ*854>0.4] （9） 式中：[ρ*i]表示中心波长为i（单位：nm）的波段的表观反射率。 4.2 基于辐射传输方程的大气校正参数计算 Kaufman等人经过大量研究发现，由于2.1 μm波长比大部分气溶胶微粒的直径要大，故该波段受气溶胶影响可以忽略，当地表反射率较低时，即为暗像元时，2.1 μm表观反射率与0.47 μm及0.66 μm的地表反射率存在很好的线性关系： [ρ470=ρ*2 1004≈ρ2 1004ρ660=ρ*2 1002≈ρ2 1002] （10） 式中：[ρ*2 100]表示中心波长为2 100 nm的波段的表观反射率；[ρ470，ρ660]表示中心波长分别为470 nm，660 nm处的地表反射率。因此根据式（10）可以得出暗像元在470 nm和660 nm处的地物反射率。在MODTRAN模型中，在校正参数一定的情况下，反射率与辐亮度有如下关系： [L=Lp+（Gb+Gt）ρ1-Sρ] （11） 式中：[L]表示某一波段传感器接收到的辐亮度；[Lp]表示对应波段的程辐射亮度；[S]为大气球面反照率；[ρ]为地表反射率，[Gb]，[Gt]为辐射校正参数，[Gb+Gt]可以通过能见度查找表查找出[Lp]，[S]，根据式（11）以及在470 nm和660 nm处的地物反射率和辐亮度来计算。 5 大气辐射校正 当考虑邻近效应且地面为朗伯面时，反射率与辐亮度有如下关系： [L=Lp+Q1-SQeGb+Qe1-SQeGt] （12） 式中：[Lp]表示对应波段的程辐射亮度；[S]为大气球面反照率；[Q]为目标反射率；[Qe]为背景反射率。 利用水汽含量反演模块和能见度反演模块中的方法的到的水汽含量和能见度建立反射率查找表，根据式（12）的关系完成反射率与辐亮度的转换。 大气辐射校正过程如下： （1） 将已得到的平均水汽含量、能见度作为已知建立反射率查找表； （2） 通过查找表求解出式（8）中的参数：[Lp]，[S]，[Gt]，[Gb]等。 （3） 用计算到的参数和式（8）校正高光谱图像得到大气辐射校正后结果图像。 6 实验分析 为了验证所设计的大气校正模块的准确性，将在东天山土屋—延东地区拍摄的高光谱数据作为实验对象，做大气辐射校正。把处理后结果与地面拍摄的实际的光谱曲线和其他软件的校正结果相比较，计算同一坐标下各校正结果的光谱特征，验证本模块的有效性。 本次遥感数据获取使用了澳大利亚HyVista公司研制的HyMap成像光谱仪，于2002年10月10日格林威治时间8时28分30秒获取东天山土屋?延东地区（北纬42.156 174°，东经94.245 092°，平均地面海拔910.433 960 m，飞行高度为3 377.978 153 m）的航空成像光谱信息，波长范围为可见光、近红外（VNIR）、短波红外（SWIR）和中热红外（MIR和TIR），波段共124个，光谱分辨率约15 nm，图像大小651×11 866，像元大小4.9 m×4.9 m，原始数据大小约1.3 GB。该遥感数据还提供了地物真实反射率数据，为评价大气辐射校正结果提供了有效的参考。 6.1 伪彩色合成对比 利用本方法校正后的图像与反射率图像以及ENVI中的FLAASH模块校正后的图像对比如下： 图4 校正得到的伪彩色图像对比结果 6.2 对应坐标光谱特征比较 校正结果与ENVI中FLAASH模块校正结果比较如图5所示。校正结果与地面实际光谱比较如图6所示。 在校正前后的结果图像中，同一坐标下两个光谱曲线的光谱角和光谱相关系数如表1所示。在表1中“A”表示ENVI的处理结果与地面测量出的辐亮度光谱曲线做比较，“B”表示本文所用方法的处理结果与地面测量出的辐亮度光谱曲线做比较。从表1可以看出，本项目校正后的光谱与地面实际光谱相比，光谱角在0.1～0.24之间，相关系数在0.90～0.99之间。说明校正结果稳定，误差在可接受范围内，与地面实际光谱较相识，达到了预期目标。 表1 校正结果与地面反射率光谱比较 图5 对应图像坐标处ENVI校正结果与本系统校正结果比较 同时，我们注意到校正过程中存在难免的误差，为此做如下可能性分析：系统如果在调用MODTRAN时建立的查找表数目偏少，会造成一定的误差，如水汽反演过程中得到了经过第一次迭代计算的水汽含量，需要在水汽查找表中查找对应的程辐射，通过线性插值得到，而如果建立更加细致的水汽含量查找表则减小误差，但必须付出更多的时间代价；另外在水汽含量迭代计算过程开始时有一个水汽含量初值的选择，若初值选择不当则会使误差扩大；在能见度反演模块中，将水汽含量的平均值作为每一个像素点下的水汽含量，而在一幅高光谱图像中，如果水汽含量分布极不均匀，那么必定会造成一定的误差，因此在必要时，需要将图像拆分成更小的数据块进行单独的大气校正。 图6 对应图像坐标处本系统校正后光谱和地面拍摄反射率光谱 7 结 语 从以上实验可以得出本系统的结果与ENVI校正结果较为吻合，校正后辐亮度值与ENVI的结果或者地面拍摄数据相比都比较吻合，只有在受水汽影响的波段误差偏大，在实际应用时只需要去除受水汽影响严重的波段即可。而且利用并行加速，可以大大减少辐射校正的运行时间，提高运行效率，使大气辐射校正达到了既快速又精确的目标。最后通过光谱特征分析，验证了大气辐射校正的有效性。 参考文献 [1] JENSEN J R.遥感数字影像处理导论[M].陈晓玲，译.北京：机械工业出版社，2007. [2] FRASER R S， FERRARE R A， KAUFMAN Y J， et al. Algorithm for atmospheric corrections of aircraft and satellite imagery [J]. Int J Remote Sens， 1992， 13： 541?557. [3] 亓雪勇，田庆久.光学遥感大气校正研究进展[J].国土资源遥感，2005（4）：1?6. [4] 丁丽霞，周斌，王人潮.遥感监测中5种相对辐射校正方法研究[J].浙江大学学报：农业与生命科学版，2005，31（3）：269?276. [5] 陶东兴，赵慧洁，贾国瑞.基于MODTRAN的双查找表法反演高光谱MODTRAN数据的水汽含量[J].光电学报，2012，32（10）：1?6. [6] 赵祥，梁顺林，刘素红，等.高光谱遥感数据的改正暗目标大气校正方法研究[J].地球科学，2007，37（12）：1653?1659. [7] 彭妮娜，易维宁，方勇华.400～1 000 nm波段反演气溶胶光学厚度的暗像元法[J].红外与激光工程，2008，37（5）：878?883. [8] SCHLAPFER D， BOREL C C， KELLER J， et al. Atmospheric precorrected differential absorption technique to retrieve colummar water vapor [J]. Remote Sens Environ，1998，65： 353?366. [9] KAUFMAN Y J， SENDRA C. Algorithm for automatic atmospheric corrections to visible and near?IR satellite imagery [J]. Int J Remote Sens，1988， 9（8）： 1357?1381. [10] KAUFMAN Y J. Atmospherically resistant vegetation index （ARVI） for EOS?MODIS [J]. IEEE Transactions on Geoscience and Remote Sensing ，1992， 30（2）： 261?270. |
| 随便看 |
|
科学优质学术资源、百科知识分享平台，免费提供知识科普、生活经验分享、中外学术论文、各类范文、学术文献、教学资料、学术期刊、会议、报纸、杂志、工具书等各类资源检索、在线阅读和软件app下载服务。