使用MODIS Level 1B 1KM 数据反演AOD实验流程

1、实验思路

实验的整体思路：
第一部分是MODIS数据预处理工作。

数据预处理工作主要包括对 MODIS LIB数据进行去除蝴蝶结效应、几何校正、波段合成、裁剪等操作

第二部分为城市地区AOD及PM2.5质量浓度反演部分。
（1）AOD反演

云监测：首先对 MODIS LIB数据进行云检测，剔除云层较厚的地区；
地表反射率的确定：然后进行暗像元识别，以及确定暗像元红蓝波段与短波红外波段地表反射率的关系
查找表的构建：同时通过6S辐射传输模型生成适用于城市地区的查找表，找表文件是IDL程序调动6S模型生成的；
基于以上工作进行AOD反演

（2）在此基础上进行PM2.5质量浓度反演

第三部分为反演结果的分析验证与精度评价。

通过结合地面监测站点的实时监测数据，对反演结果在时间、空间、强度方面进行分析，并对反演结果进行精度验证及评价

暗像元法：

2.1 μm 的短波红外波段由于受大气气溶胶的影响较小，可以忽略不计，使得其表观反射率和地表反射率近似相等。
Kaufman 等人通过大量的飞机观测实验发现红、蓝和短波红外通道的地表反射率具有很好的线性关系。
因此，可以从短波红外波段的表观反射率获得红蓝波段的地表反射率，然后从红蓝波段的表观反射率去除地表贡献，获得大气参数 S、ρ0、T( μs) T( μv) ，进而得到 AOD。

2、实验数据

主要使用到的MOD021Km HDF数据集是发射率、反射率、角度

2.1 数据介绍

依据标准数据产品分级系统的规定， MODIS数据又可按处理情况划分为五个等级，即0、1、2、3、4级数据。

0级数据：即 MODIS原始影像资料，是指卫星地面接收站获取到影像之后，通过格式化同步、去除干扰、遥感纠错以及变换格式等处理步骤，得到的包括全部数据信息的原始数据，其扩展名称为（PDS）。
1级数据：即已经被赋予了定标参数的数据。对原始影像资料进行处理，借助辅助数据资料进行注释解译，并进一步经过时间配准等计算和增补步骤之后得到1级数据。1级数据主要存储为HDF格式，即数据文件扩展名为（HDF）。
2级数据：2级数据是将1级数据进一步处理加工后得到的，与1级数据的覆盖范围一致而且拥有同样的空间分辨率
3级数据：3级数据是在2级数据的基础上，对传感器成像过程中产生的蝴蝶结效应做了处理，3级数据一般情况下具备一致的时空描述变量并且拥有一定的一致性与完整性，在3级数据的基础上就可以集中进行来自单一技术的观测方程以及通用模型等方面的研究工作了。
4级数据：4级数据是根据参数文件所提供的参数，对遥感影像几何纠正、辐射校正，从而令影像的每一个像素都拥有准确的地理信息、反射率信息、辐射亮度信息等，即综合分析3级及以下数据而得到的应用级产品不可或缺的基础。

问题：

反演AOD为什么要用到发射率和反射率的合成数据？
反演AOD为什么要用到角度数据集？反演公式的参数、构建查找表的参数
反演AOD和大气校正的关系？
MOD02 1Km中的反射率数据是地表反射率还是表观反射率？表观
MOD02 1Km中的质量数据是干什么用的，体现哪些方面的质量？

2.2 数据获取

https://ladsweb.modaps.eosdis.nasa.gov/search/

MOD021KM.A2012156.0255.005.2012156094454.hdf

3、实验操作

实验操作的具体流程：

3.1 打开数据-设置实验的输入输出路径，并打开原始文件MODIS Level 1B 1Km，进行辐射校正：

使用ENVI打开MODIS HDF数据的方式不一样，能在layer manager中显示的波段也不同，详细了解可以参照：ENVI打开MODIS HDF数据的几种方式及对比。
关于MOD02 1Km数据的官方介绍是：MOD021KM - Level 1B Calibrated Radiances - 1km

1、设置在file reference directory 里面设置好输入输出路径
2、打开HDF文件，file-open as-EOS-MODIS，这个是MODIS L1B（MOD02）这个级别数据的打开
3、打开后点击zoom to full extent,在窗口显示全部的图像范围
4、打开data manager，可以看到一共打开了6个数据集

Emissive发射率数据集，
Radiance辐射亮度的数据集，
Reflectance反射率数据集，
Longitude
Latitude
Quality

其中，Emissive发射率数据集， Radiance辐射率的数据集， Reflectance反射率数据集，都进行过了辐射定标，即这种打开方式，图像灰度具有明确的物理含义，不需要再进行波段运算进行定标。虽然这种打开方式速度快，但是适用的数据有限，打开后得到的图像波段也有限。比如MOD02数据中也有太阳/传感器天顶角、方位角波段，用这种方式就无法打开。

5、打开cursor value，可以看到DN值小于1，可以看出，ENVI打开MODIS数据时，自动做了辐射校正，DN值是定标后的DN值；这个在ENVI CLASSIC里面可以看到，它调用的是ENVICLASSIC里面的算法，在ENVI CLASSIC 里面打开file- reference-Miscellaneous [ˌmɪsəˈleɪniəs]，其中的一个选项：Auto-Correct ASTER/MODIS。如果想要查看原始的DN值那么，在打开图像之前需要在ENVI CLASSAIC中把Auto-Correct ASTER/MODIS旁边的Yes改为No。
6、加载经纬度数据集，Load Grayscale,打开的经纬度数据是灰度图，它的行列号以及分辨率和影像相同，它的DN值就是该像元的经度或者纬度。
所以HDF是自带有地理定位信息的数据。

概念

发射率也叫比辐射率或发射系数，是指地物发射的辐射通量与同温度下黑体辐射通量之比。地物的发射率与地物的性质、表面状况(如粗糙度、颜色等)有关，且是温度和波长的函数。
辐射亮度简称辐亮度，辐射亮度表示面辐射源上某点在一定方向上的辐射强弱的物理量。
反射率是物体反射的辐射能量占总辐射能量的百分比。地表反射率：地面反射辐射量与入射辐射量之比，表征地面对太阳辐射的吸收和反射能力。反射率越大，地面吸收太阳辐射越少；反射率越小，地面吸收太阳辐射越多。表观反射率：表观反射率就是指大气层顶的反射率,辐射定标的结果之一，大气层顶表观反射率，简称表观反射率，又称视反射率。（=地表反射率+大气反射率。所以需要大气校正为地表反射率）。
参考：https://blog.csdn.net/lucky51222/article/details/22610771

3.2 发射率数据的几何校正

7、发射率数据的几何校正，打开右侧Tool Box目录中的Geometric Correction-Georeference by Senser，找到Georeference MODIS双击，选择发射率数据集，点击OK

然后，弹出Georeference MODIS Parameters对话框，设置要输出的坐标系，选择UTM，基准面Datum选择WGS-84；带号选择研究区域所在的带号，如果不知道带号的话，点击右侧的Set Zone，输入经纬度自动计算出带号。

然后先不要点击OK，先输出GCP文件，它会在它自带的经纬度文件中生成一些控制点，然后将输出的文件命名为GCP.pts，设置好输出路径。最后是修正双眼皮效应（蝴蝶结效应），选择yes，最后点击OK。然后出现Registration Parameters对话框中，注意分辨率是1000（默认）选择保存路径点击OK即可。

点击Ok，（这种自带定位信息的图像进行几何校正比较简单）

最后输出校正结果：

此外进行MODIS数据的几何校正还有另外一种方法：Registration-Reprojecrt GLT with Bowtie Correction, 但是这种方法不能输出GCP地面控制点

3.3反射率数据的几何校正（相同方法进行）

与发射率数据校正方法一样，也可以使用Wrap from GCPs： Image to Map Registration工具进行校正。投影，基准面等设置照旧，但是不用再设置GCP输出路径了，然后设置输出路径，点击OK

反射率数据几何校正结果：

3.4 合成与裁剪-发射率数据几何校正结果与反射率数据几何校正结果进行合成与裁剪

9、然后打开data manager，可以看到发射率和反射率几何校正的结果是多波段的数据，下面对这两个多波段文件进行波段的叠加，就是让他们合成一个多波段的文件，这样便于后面的分析。

首先找到layer stacking工具

选择反射率几何校正结果和发射率几何校正结果，然后点击空间裁剪按钮，Spatial Subset，

这里给出了四种裁剪方式
Image方式是手动地在图像上选择一个矩形的一个范围

Map是给出左上角和右下角的经纬度坐标，是基于坐标的裁剪

第三，file，是基于另外一个文件进行裁剪
第四，ROI/VF，是感兴趣区域或者是矢量裁剪；
在这里，我们选择第四种方式进行空间裁剪：选择一个外部的矢量文件，beijing.evf,这是ENVI矢量格式，下图显示的是可以用shapfile转换成ENVI的格式，选择beijing.shp，点击OK

自动确定了选择范围，这个范围是不规则矢量边框的范围
点击两次OK后会出现这个对话框：

注意这个顺序有问题，应该是反射率在上，发射率在下，怎么调换呢？点击Reorder Files，然后点击鼠标中键即可

点击OK，再设置输出路径

结果（这个工具只能对有地理坐标数据的多波段数据进行合成叠加，如果没有地理坐标的话只能用其他的工具）

3.5 打开角度数据

10、查看角度数据的属性信息。搜索View HDF Dataset Attributes工具。View HDF Dataset Attributes用于打开HDF图像的基本属性

打开HDF文件，然后选择传感器的方位、天顶角和太阳方位角和天顶角这四个角度属性，点击OK

得到只包含选定的四个属性的信息的对话框

要想打开他们，需要在ENVI菜单栏中，file-open as-generic formats-HDF4打开HDF Dataset Selection对话框，

打开四个角度文件，点击OK，这样就把四个角度信息在layer manager中打开了

3.6 角度数据的合成

在data manager中可以看到，这四个文件是单波段文件，下面要把他们合成为一个文件，让这四个文件作为四个波段在一个多波段文件中，这样便于整体处理。此处的合成不能使用layer stacking 进行合成了，使用New File Builder工具，这个工具可以对没有地理坐标的数据进行波段合成

选中四个角度文件后点击OK，这一步先不做空间裁剪，这个角度文件也要按照一定的顺序：卫星天顶角、方位角，太阳天顶角、方位角，data set顺序分别为14、15、17、18

合成结果：

3.7 角度数据重采样

11、下面是角度数据的几何校正，因为他没有地理坐标，所以进行几何校正的时候可以用上一步发射率数据得到的控制点文件对其进行几何校正，但是我们来看一下，角度数据的行列数是271-406，分辨率比较低，而原始的反射率图像的行列数是1354-2030，所以需要对角度数据集进行重采样，采样结果要和反射率数据集的分辨率/或者行列数相同

打开数据重采样工具Resize Data，选中角度合成结果，点击OK，
没有分辨率信息可以通过行列号信息进行重采样

点击Layer manager-view-remove all layers，然后打开data manager,选择角度重采样结果.dat，然后右键load Grayscale；选择发射率数据集（因为GCP文件是从它得来的），也Load Grayscale。可以看到，角度合成的数据集行列号已经和发射率数据一致了，下面用控制点校正的方法对角度数据进行几何校正。

3.8 角度数据的几何校正

双击warp from GCPs : Image to Map Registration工具，然后选择GCP文件；

选择过GCP文件后进行参数设置，注意把分辨率设置为1000，然后点击OK

然后选择待校正的文件角度重采样结果.dat

几何校正方法选择Triangulation；重采样方法选择Cubic Convolution

角度重采样数据几何校正结果

3.9band math运算与角度数据裁剪

角度数据的DN值现在是一个4位数字，整数。我们在View HDF Dataset Attributes工具中查看一下原始图像的角度属性，选择其中一个角度查看就可以了

可以看到，单位是度，范围是0到18000，可以猜测到角度放大了100倍，然后看到缩放系数是0.01，也就是说，角度数据集的DN值要乘以0.01，才是真实值。
搜索一下band math工具

然后直接OK

点击Map Variable to Input File

在弹出的Band math Input File对话框中选择几何校正后的角度合成数据集，然后点击OK

此时再进行空间裁剪，裁剪后的结果：

3.10云检测-气溶胶反演

将MODIS云检测和气溶胶反演的IDL拓展程序复制粘贴到ENVI安装目录下的ENVI53-extensions中，然后重启ENVI，可以看到在右侧的Tool Box中的Extensions里多了新添加的这两个工具。

首先进行云检测，点击Open打开反射率和发射率合成的结果，然后点击OK，对反射率发射率合成的数据集做云检测，去除了云的像元
反射率在上发射率在下这与算法的设置有关

云检测结果

气溶胶反演
打开Extensions里添加的气溶胶反演工具，按照对话框的提示打开不同的文件
首选打开云检测结果：

点击OK，然后打开角度数据
点击OK，然后选择查好表文件

这个查找表文件是IDL程序调动6S模型生成的
点击OK，然后设置气溶胶反演结果的保存路径

气溶胶反演结果：

然后这个图像像元的DN值就是该像元处的气溶胶光学厚度，AOD值越大，空气质量越差

3.11 后处理-重新划分颜色等级

在layer manager中选择气溶胶反演结果图层，右键，选择New Raster Color Slice，然后选择要进行密度分割的图层，点击OK

默认划分了16个等级，先全部删除，然后手动设置他的等级个数以及对应的颜色，我们从外部导入密度分割的文件，“气溶胶颜色等级.dsr”将其分为7个等级，然后点击OK，得到重新分级后的结果

然后点击layer manager中的Slices，右键，export color slices – class image

设置完输出路径点击Ok。

3.12 后处理-做研究区域的掩膜

先打开北京市的矢量文件beijing.shp

然后选择Apply Mask

输入文件选择气溶胶颜色分级显示的结果，然后再select mask band 这里选择mask options-build masks，然后弹出下面这个对话框

它自动有一个0到0的等级，我们把他删掉，然后options下面选择import EVFs，然后选择beijing.shp，点击OK

在弹出的对话框中选择栅格文件来确定掩膜文件的行列号和分辨率。掩膜文件的应用有一个前提，你要将其应用于某个影像，就要和该影像的分辨率和行列号保持一致，所以，建议掩膜文件的时候，就选择他要应用于的那个影像文件

点击OK后出现下面的对话框，并在设置好输出路径和文件名后点击Ok

掩膜建好之后就选择到mask band上来，之后点击OK，然后设置掩膜参数，mask value 设置为-1，就是不要用有效值，因为AOD一般是一个大于0的值，所以可以设置为-1

很不巧，背景处理后，-1的区域变成了红色，那我们使用 edit envi header

得到最终的结果

参考文章：
http://blog.sina.com.cn/s/blog_764b1e9d01015oq5.html
http://muchong.com/t-10968593-1-authorid-3088580
https://www.jianshu.com/p/559738c230a2
https://wenku.baidu.com/view/b45c09f09e31433239689336.html