操作-大气校正,辐射定标,气溶胶反演

基于RS\GIS监测洪灾变化上机操作实例

基本原理：

①大气校正

遥感图像在获取过程中，受到大气吸收与散射、传感器定标、地形等因素的影响，且会随时间的不同而有所差异。利用多时相遥感图像的光谱信息检测地物变化的重要前提是要消除不变地物的辐射值差异。

大气校正的目的是消除大气和光照等因素对地物反射的影响，大多数情况下，大气校正是反演地物真实反射率的过程。

目前可以进行大气校正的模块有很多种，如最早的MODTRAN 4+，6S (Second Simulation of the Satellite Signal in the Solar Spectrum)，ACORN，ATREM，在ERDAS IMAGINE 8.7上的模块ATCOR，以及ENVI上的模块FLAASH（基于MODTRAN）。

FLAASH可对LANDSAT，SPOT，A VHRR，ASTER，MODIS，MERIS，AATSR，IRS等多光谱、高光谱数据、航空影像及自定义格式的高光谱影像进行快速大气校正分析。下面的大气纠正步骤，都是基于FLAASH进行的。

②辐射定标

当我们拿到一幅原始影像，先要进行辐射定标，目的是把图像上的DN（Digital Number）值转为辐亮度或者是反射率。辐射定标的结果可以是表观辐亮度（L），也可以是表观反射率（ρ）。

计算表观辐亮度（L）的公式为：

Radiance=（（Lmax-Lmin）/（Qcalmax-Qcalmin）*（Qcal-Qcalmin）+Lmin ①

其中：Radiance 是表观辐亮度，注意单位是W/m2·sr·μm；Qcal为像元DN 值（也就是影像数据本身）；Qcalmax为传感器处最大辐亮度值所对应的DN值，一般为255；Qcalmin 为传感器处最大辐亮度值所对应的DN值，一般为0；Lmax 和Lmin是从参数表中查询，Lmin为光谱辐亮度的最小值，单位同L；Lmax为光谱辐亮度的最大值，单位同L。

计算表观反射率（ρ）的公式为：

ρ =π*L*d2/（ESUN*cos（θ））②

其中：ρ为表观反射率；L为①式中计算出来的表观辐亮度；d为日地距离；ESUN为大气层外的太阳辐射，也可以说是传感器接收处的太阳辐射；θ为太阳天顶角（这个可以通过影像的元数据获取）。以上参数可以查询下表获得。

对于TM影像来说，当LANDSAT传感器获取影像以后（Level 0），会将其转化为32 位浮点型的绝对辐亮度。之后进一步处理，将绝对辐亮度变为8位的DN 值，这也就是我们购买后拿到的数据（Level 1）。

如果要将L1的DN值转化为传感器处的辐亮度值（at-sensor spectral radiance），对于TM是8bit来说，Qcalmax取255，Qcalmin取0，可以将公式①简化为下面这个公式：

上面的这个公式也可以改为：

其中，

Grecale即为gain，Brescale即为offset。各个波段的，以及

和见下表。

需要注意的是，上述参数在2003年5月5日前后是不一致的，所以在操作时，一定要搞清楚影像获取的时间。ETM+影像与TM影像的Lmin、Lmax、Grescale、Brescale、d及ESUN的参数不一样，可以查询下列表格：

③气溶胶反演

我们采用FLAASH模块中的暗目标法来反演气溶胶光学厚度。暗目标法是由Kaufmaun提出，它是利用660nm和2100nm处的反射率估算气溶胶量，由于2100nm 波长比大部分气溶胶微粒的直径要大，故该波段受气溶胶影响可以忽略；在大量的试验中，他发现2100nm的植被反射率与660nm植被反射率之间存在稳定的关系，因此可以直接利用2100nm的植被反射率来获取660nm处的植被反射率。气溶胶的影响会使得实际获得的植被反射率与理论反射率存在一定差异，FLAASH模块中正是利用了这个差异来反演气溶胶的光学厚度值。

要在FLAASH模块中获取图像气溶胶含量，传感器必须具有660nm和2100nm 附近的通道，这些通道主要是用于获取“黑暗像元”，条件为，如果输入图像中还具有800nm 和420nm 附近的通道，可以用于消除阴影和水体，条件为。

数据准备：

经过几何粗校正、没有经过大气校正的Landsat-TM5（30米分辨率）遥感影像两景，分别位于2010-7-2与2010-9-4文件夹中。

其中有关影像的相关信息存放于后缀为_MTL.txt的文件中，TM影像1-5及7波段的中心波长信息，存放于中心波长.txt文件中，用于大气校正。

处理过程：

①大气校正：打开ENVI软件，点击主菜单栏的file下的open image file，打开2010-9-4中的LT51520412010247KHC00_B1—B5及B7影像，再选择主菜单栏的Basic Tools---Preprocessing---Calibration Utilities---Landsat Calibration，对6个波段进行辐射定标。

在Landsat Calibration Input File列表中，先选择LT51520412010247KHC00_B1第一波段影像进行处理，点击OK。

参考处理影像相关信息文件LT51520412010247KHC00_MTL.txt修改ENVI Landsat Calibration对话框参数，如下图所示，输出文件命名为band1_ca.TIF。按照相同方法依次对其余5个波段进行辐射定标，处进行波段选择时需要注意。

选择主菜单栏的File---Save File As---ENVI Standard，点击Import Files，导入辐射定标后的band1_ca.TIF—band5_ca.TIF及band7_ca.TIF 6个波段文件，再点击Reorder Files对6个波段从1至7依次进行排序，最终保存为的20100904_ca文件包含6个波段。

接下来进行FLAASH校正操作。选择主菜单栏的Spectral---FLAASH，出现FLAASH Atmospheric Correction Model Input Parameters窗口，FLAASH 模块的操作界面分为三块：最上部设定输入输出文件；中间设定传感器的参数；下部设定大气参数。

首先设定输入输出文件。FLAASH 模块要求输入辐亮度图像，输出反射率图像。之前我们进行了辐射定标，得到辐亮度图像，在这里要把BSQ 格式的图像转换为BIL 或者BIP 格式的图像，选择主菜单栏的Basic Tools---Convert Data(BSQ，BIL，BIP)。然后再Input Radiance Image 中选择转换格式后的图像。

当输入图像后，需要导入已经准备好中心波长.txt文件，其中含有一列TM 每个波段中心波长的信息。然后需要选择Scale Factor，即原始辐亮度单位与ENVI 默认辐亮度单位之间的比例。ENVI 默认的辐亮度单位是μW/cm2·sr·nm，而之前我们做辐射定标时单位是W/m2·sr·μm，二者之间转换的比例是10，因此在下图中选择Single scale factor，填写10.000。

在Output Reflectance File里面设定输出文件的文件名和位置。

设定传感器参数。首先是Scene Center Location，即遥感图像中心的坐标，以及Flight Date、Flight Time GMT（the Greenwich Mean Time），后两者可以在TM 的LT51520412010247KHC00_MTL.txt文件中找到，填入即可。遥感图像中心坐标可以在TM的任一波段的头文件中获取到，在任一波段文件处右键点击Edit Header，进入后点击Edit Attributes---Map Info，点击，可以选择DMS或DDEG中任一种形式填写入FLAASH窗口。

在Sensor Type 菜单中选择Landsat TM5，此时Sensor altitude 自动填上为705km。而Pixel Size 填为30m。

根据遥感影像研究区实际情况，填写Ground Elevation，因为巴基斯坦南部地区平均海拔在20m左右，所以填写为0.02km。

最关键的为大气参数部分:

a)Atmospheric Model(大气模式): 共有Sub-Arctic Winter (SAW)，Mid-Latitude

Winter (MLW)，U.S. Standard (US)，Sub-Arctic Summer(SAS)，Mid-Latitude Summer (MLS) 和Tropical (T)。根据经纬度和时间可以选定研究区的大气模式，见下表，这里选Tropical (T)。（研究属于北纬28°）

b)Aerosol Model（气溶胶模式）：有Rural，Urban，Maritime和Tropospheric四种

选择。根据实际情况选择即可。关于此四种模式的解释见下图。

c)这里选Maritime。

c) 当我们选择TM 时，可选的参数还有Aerosol Retrieval 和Initial Visibility。这两

个参数对最后的结果有相当重要的影响，因此最好能调查到当地的Initial Visibility。此外，AERONET 在全世界各地有测定AOD（Atmospheric Optical Depth）的站点，可以查询AOD 以后转换为消光系数，通过消光系数估算能见度，此步骤比较繁琐，在此不予详述。如果采用Aerosol Retrieval 中的K-T算法计算Visibility，且能够计算出结果的话，则采用K-T 算法的能见度。

d) 关于Aerosol Retrieval。如果选择了下拉菜单中的K-T method，那么需要在Multispectral Settings 中设定参数，在Kaufman-Tanre Aerosol Retrieval/Assign Default Values Based on Retrieval Conditions 中选择Over-land Retrieval Standard (660:2100nm)即可。根据不同的研究区可以设定不同的模式。

这些参数用于确定黑暗像元，用于气溶胶反演；

KT Upper Channel：建议选择2100nm附近的通道；

KT Lower Channel：建议选择660nm附近的通道；

Maximum Upper Channel Reflectance：建议设置为0.1，即：。

Reflectance Ratio：为反射率比值，建议设置为0.45，即：

Cirrus Channel (optional)：确定云的通道，建议设置为1367—1383nm左右的通道。

水汽、气溶胶和云反演通道设置波长范围如下图所示：

上面这些参数波长范围的选择可以参考中心波长.txt文件，设置完Apply即可。

水量结果。

与裸地波谱曲线更能真实地反映情况。

对2010-9-4影像进行大气校正后，根据相同原理及步骤对2010-7-2影像进行大气校正（叙述省略，可以参照截图）。

②提取水体范围：运用下式提取巴基斯坦洪灾部分区域的水体面积。

MNDWI=（TM2- TM5）/（TM2+ TM5）

选择主菜单栏的Basic Tools---Band Math，点击Import Files，输入计算公式：

float((b1-b2))/float((b1+b2))

点击OK 。其中b1选择TM2波段，b2选择TM5波段，命名保存至文件路径，点击OK 进行计算。

通过计算得到的影像中，水体为正值，在0-1之间，土壤、植被及建筑物等物体为负值。

接下来继续使用Band Math功能对提取MNDWI后的影像进行二值化，提取出水体范围。因为水体的值大于等于0，而其他地物的值小于0，所以可以使用下列语句对水体赋值为1，而其他地物赋值为0：

(b1 ge 0)*1+(b1 lt 0)*0

其中ge、lt分别表示“大于等于”和“小于”。

点击OK。其中b1选择计算好的MNDWI波段，命名保存至文件路径，点击OK进行计算。

计算得到二值化后的影像。

按照相同原理和方法对2010-7-2影像进行MNDWI提取和二值化，步骤省略。

因为影像最左侧边缘有影像质量问题引起的干扰，所以需要对其进行裁剪，以免造成误差。这里我们选用感兴趣区进行裁剪的方法。选择主菜单栏的Basic

Tools---Region Of Interest---ROI Tool。

我们可以较为粗略的画一个框覆盖住有误差的地区，画完后在框内点击右键结