1 前期准备与本文理论部分

1.1 几句闲谈

  前面几篇博客介绍了基于Landsat这一多光谱遥感图像数据的多种地表温度（LST）反演方法，大家可以参考博客1、博客2、博客3；那么接下来，我们就将基于比多光谱数据可以说是更进一步的高光谱卫星数据——大名鼎鼎的Hyperion数据，进行多种其他地表参数的反演。其中，在此之前有必要先了解一下国内外主流的星载高光谱传感器及其平台的相关信息。
  好啦，了解相关背景之后，话不多说，我们开始本次的内容。

  其中，本文所用遥感影像相关数据：
  点击链接：https://pan.baidu.com/s/16JUtCFhXlprRLjIhNBlikA
  提取码：csdn

1.2 背景知识

1.2.1 Hyperion数据介绍

  本文所使用的遥感影像数据不同于上述三篇博客所用的Landsat-7 ETM+影像数据，其来自地球观测卫星-1（Earth Observing-1，EO-1）。EO-1是美国国家航空航天局（National Aeronautics and Space Administration，NASA）新千年计划（New Millennium Program，NMP）地球探测部分中第一颗对地观测卫星，其目的即为在21世纪接替Landsat-7卫星，于2000年11月发射升空。除此之外，NMP目前还包括深空探测（Deep Space，DS）、空间技术（Space Technology，ST）两个部分。
  EO-1卫星轨道参数与Landsat-7较为近似，以期实现两颗卫星图像每天具有1至4景的重叠，从而进行二者的对比。此外，EO-1已于2017年2月停止服役。
  EO-1搭载了三种传感器，分别为高光谱成像光谱仪（Hyperion）、高级陆地成像仪（Advanced Land Imager，ALI）与线性标准成像光谱仪阵列大气校正器（the Linear Etalon Imaging Spectrometer Array Atmospheric Corrector，LAC）。一般地，传统的陆地资源卫星只能提供为数不多的多光谱波段，并不能很好满足日常实际研究、运用的需要；而借助具有242个波段、光谱范围为356至2578纳米的EO-1 Hyperion传感器，可获得更具价值的高光谱数据。
  EO-1遥感影像命名格式如下：
  EO1SPPPRRRYYYYDDDXXXML_GGG_VV
  其中，EO1为EO-1卫星代号，S为所用传感器代号（H为Hyperion传感器，A为ALI传感器），PPP为成像时目标所处全球参考系统（Worldwide Reference System，WRS）的轨道（Path），RRR为成像时目标所处WRS的行（Row），YYYY为成像年份，DDD为成像当日在该年份的天数，XXX分别为Hyperion、ALI与AC三种传感器的开关状态（1为开启，0为关闭），M为指向模式【（N为天底模式（Nadir），P为点在轨道模式（Pointed Within Path/Row），K为点不在轨道模式（Pointed Outside Path/Row）】，L为图像长度【F为全景（Full Scene），P为局部景（Partial Scene），Q为次级局部景（Second Partial Scene），S为样例（Swatch）】，GGG为影像地面接收站，VV为影像版本编号。
  一般地，遥感卫星传感器主要有两大类型：摆扫式（Whisk Broom Scanners）与推扫式（Push Broom Scanners）；Hyperion属于后者。其242个波段分为可见光近红外波段（V-NIR）与短波红外波段（SWIR）；其中，1至70波段属于V-NIR通道，71至242波段属于SWIR通道。两个波段之间具有20个波段的波长数值相互重叠，其分别用两套不同的敏感元件收集各自信号。
Hyperion产品波段信息如下所示。

  一般地，辐射校正包括辐射定标、大气校正和太阳及地形校正，用来消除辐射误差；而上述“辐射校正”包括正射校正，即使用地形数据的几何校正，不包括大气校正。
  Hyperion产品分为两级，Level 0与Level 1。前者为原始数据，其仅用来生产Level 1产品。Level 1产品则可以继续分为L1A、L1B、L1R、L1Gs与L1Gst等。其中，L1B产品与L1R产品分别由美国TRW与USGS处理生成。上述两种产品与L1A产品的最大不同在于，前二者纠正了V-NIR波段与SWIR波段的空间错位问题。
  而通过图像原始数据中“README.txt”文件可知，本文所使用的数据为L1Gst级别（大家自己操作时，也依据自己文件的实际情况判断即可）；如以上内容所述，其已经过部分预处理过程，如正射校正、坐标转换等。
  Hyperion产品图像数据的空间分辨率为30米。
  如上述内容所提到的，Hyperion产品共有242个波段，波长覆盖范围为356至2578纳米。而在这其中，由于紫色光与短波红外两端有些通道的响应度较低，数据利用价值不大，因此系统并未对这些波段加以辐射校正——即对于正式分发的Hyperion产品Level 1数据，上述242个波段中，1至7波段、58至76波段、225至242波段亦称之坏波段（Bad Bands），其数值均为0值。除此之外，121至126（或127）波段、167至178（或180）波段、224（或222至224）波段由于受到水汽吸收影响较为严重，亦需要作为坏波段处理。最后，通过上述筛选后剩余的波段中，还包括波长重叠的56至57波段和77至78波段；考虑到77至78波段噪声相对前者较大，信噪比较低，因此将其同样剔除。
  上述波段筛选结果如下所示。

  此外，除上述因未定标和水汽吸收影响被剔除的波段外，还包括坏线、条纹与“Smile”效应等非正常像元。其中，坏线是指一行或一列DN值为零或非常小的像元；条纹是指像元DN值不为零但较小，与周围有明显差异的带状现象；“Smile”效应是指由于前期光谱定标而产生的光谱差异。若对定量遥感结果精度要求较高，上述这些非正常像元或现象均需要在预处理阶段加以解决。
  目前，针对Hyperion产品图像的坏线处理主要通过均值邻域法，即用坏线周围的像元像素值求平均后替换坏线中的像素——这依赖于坏线像元与周围像元的高度相关关系；若当坏线出现在地形比较复杂的区域，采用这种方法则会导致处理结果精度较低。
  对此，有国内学者提出了一种基于地物类型的坏线修复方法。该方法假设各波段图像成像在一瞬间完成，或成像过程中天气未发生变化，则不同波段图像中同一类地物的DN值应当相同。当某一波段的图像A的地物M出现部分坏线时，搜索另一对应位置无坏线的波段对应的图像B，确定B中M位置对应的DN值α；随后继续在B中搜索，找到所有DN值为α的像元N。随后返回A，求取A中N像元对应位置的DN值，将其求平均后作为M位置处的DN值。这一坏线修复方法大大缩小了修复结果与真实值之间的均方根误差，效果较好。
  条纹尤多出现于SWIR波段，其严重影响图像质量；因此，条纹去除同样十分重要。针对这一问题，有学者提出了一种全局去除条纹的方法。其原理公式如下：

  其中，〖x^’〗_ijk为第k波段对应图像第i行，第j列像素校正后的数值，x_ijk为其原值；s_ik为第k波段对应图像第i列标准差，m_ik为第k波段对应图像第i列平均值，(s_ik ) ̅与(m_ik ) ̅分别为第k波段对应图像第i列标准差与平均值的对应平均值。实际计算时，用整幅图像的平均值与标准差代替每一列的平均值与标准差。
  针对条纹去除结果的检验，可以利用最低噪声分离旋转（Minimum Noise Fraction Rotation，MNF）转换加以实现。
  “Smile”效应是指在垂直飞行方向上，像元的波长随着图像中心位置向两旁移动，产生的偏移。其中，V-NIR波段的“Smile”效应较为明显，约在2.6至3.6纳米范围内；SWIR波段的这一效应表现则不如前者明显，其约在0.40至0.97纳米范围内。“Smile”效应同样可以借助光谱信息的MNF转换加以检测。“Smile”效应往往会对植被遥感产生影响；可以通过大气纠正来降低“Smile”效应的影响。
  Hyperion产品Level 1数据以有符号的整型数据格式记录其信息，数值范围为-32768至32767；而实际地物辐射数值都比较小，如V-NIR波段最大辐射约为750 W/m2/sr/Lm，而SWIR波段最大辐射仅有350 W/m2/sr/Lm。因此，在产品生成时，系统采用了扩大因子，即分别将V-NIR与SWIR所对应波段的数值扩大为40和80倍。出于这个原因，在应用Hyperion产品Level 1数据时，需要将像元值转换为绝对辐射值。其中，两通道对应波段的计算公式有所不同：

  其中，L_ℷVNIR为V-NIR波段定标后的绝对辐射值，L_ℷSWIR为SWIR波段定标后的绝对辐射值。
  此外，在计算时需要注意，由于经过波段筛选后的图像波段不再完全连续，会出现一些间段区域。因此需要首先将同一通道内波段编号相连的波段合并，再将同一通道内全部波段合成；第二次合成后，将两个通道对应的波段分别带入上述两个公式，计算结束后再将两个通道的176个波段最终合成为一幅图像。
  编辑头文件时，除需注意中心波长的修正外，还要注意半高宽（又称半峰全宽，Full Width at Half Maxima，FWHM）的输入。若不进行大气校正或进行大气校正但不采取FLAASH大气校正方法，则无需导入FWHM这一参数。
  需要注意的是，本文重点在于Hyperion的反演操作，因此上述Hyperion图像的坏线修复、条纹去除等处理暂未写入本文。具体操作我将在后期的博客中展现。

1.2.2 遥感图像分类方法

  本文使用监督分类（又称训练分类，Supervised Classification）方法提取影像中各个不同地物部分，尤其是太湖湖面水体区域。而和前期博客一致，本文的监督分类部分依然借助ERDAS软件，选择采取平行六面体规则（Parallelepied）这一非参数规则（Non-parametric Rule）与最大似然规则（Maximum Likelihood）这一参数规则（Parametric Rule）加以实现；并且尝试了特征空间规则（Feature Space）的效果。
  监督分类是指在已掌握有足够先验知识（亦即训练场地）的情况下，根据已有训练场地提供的已知属性样本选择特征参数，并训练、建立得到对应判别函数；随后进而将图像未知类别部分的像素的值代入建立得到的判别函数，依据所选择的不同判别准则，对该样本所属的地物类别进行自动分类处理；即监督分类是一种利用已知地物属性、信息，进而对未知属性地物进行分类的方法。
  正如上述内容所提到的，常用的监督分类判别规则包括非参数规则与参数规则两个部分。其中，非参数规则包括特征空间规则与平行六面体规则；参数规则则包括最小距离规则（Minimum Distance）、马氏距离规则（Mahalanobis Distance）、最大似然规则以及波普角规则（Spectral Angle Mapper）等多种。
  ERDAS IMAGINE 2015软件中“Supervised Classification”模块可以选择的规则十分齐全。本次我们依然运用这一模块，并选择采取平行六面体规则与最大似然规则加以实现。平行六面体规则是指：根据训练样本的图像亮度值，形成一个N维的平行六面体数据空间；其余像元的光谱数值如果落在平行六面体任何一个训练样本所对应的区域，其就被划分至这一对应的类别中。其中，平行六面体的尺度是由标准差阈值所确定的，而该标准差阈值则是根据所选类的均值求出。最大似然规则是指：假设每一个波段的每一类统计都呈正态分布，并计算给定像元属于某一训练样本的似然度；随后，像元将最终被归并到似然度最大的一类样本当中。
  另一方面，遥感影像处理也可使用非监督分类（又称聚类分析或点群分类，Unsupervised Classification）加以完成。非监督分类方法是在多光谱图像中搜寻、定义其自然相似光谱集群的过程，其不必针对影像地物获取先验知识，仅依靠影像中不同类地物的光谱（或纹理）信息进行特征提取，再统计特征的差别，从而达到分类的目的；最后对已分出各个类别的实际属性进行确认。

1.2.3 大气校正

  本次实验可以不进行大气校正；但不经过大气校正，由于未消除大气影响，会使得得到的结果会有一定误差。本实验中，我们首先不进行大气校正，对叶绿素加以反演；随后再先后尝试FLAASH大气校正与QUAC快速大气校正两种方法，对得到的高光谱遥感影像加以处理，从而对叶绿素反演精度加以提升。
  大气层是介于卫星传感器与地球表层之间的一层由多种气体及气溶胶组成的介质层，其由下至上可依次分为：对流层（0至8-18千米）、平流层（8-18至50-55千米）、中间层（50-55至80-85千米）、电离层（又称热层，80-85至800千米）、散逸层（800至3000千米）。在太阳辐射到达地表再到达卫星传感器这一过程中，辐射两次经过大气，故大气对太阳辐射的作用影响较大。
  大气校正的目的即消除大气和光照等因素对地物反射的影响，广义上讲是获得地物反射率、辐射率或地表亮度温度等真实物理模型参数，狭义上讲是获取地物真实反射率数据。大气校正可以用来消除大气中水蒸气、氧气、二氧化碳、甲烷和臭氧等物质对地物反射的影响，也可以消除大气分子和气溶胶散射的影响。在大多数情况下，大气校正也是反演地物真实反射率的过程。
  有时可完全忽略遥感数据的大气影响。对某些地物分类和变化检测过程（如用最大似然法对单时相遥感数据进行分类）而言，只要影像中用于分类的训练数据具有相对一致的尺度（即均经过或未经过大气校正），大气校正并不是必需的；而当训练数据需要进行时空拓展时，其影像分类和各种变化检测则需进行大气校正。
  而对于涉及定量遥感的应用方面，如前期博客所进行的地表真实温度反演运算，大气校正则多是必不可少的环节。有研究指出，是否进行大气校正过程，对植被稀少或植被被破坏地区的NDVI计算误差影响可达50%。
  大气校正方法整体上可以分为两类：一是绝对大气校正方法，即将遥感图像的DN值转换为地表反射率、地表辐射率、地表温度等数值的方法；二是相对大气校正方法，即将原始的DN值校正为不考虑地物实际反射率情况下，相同DN值代表相同地表反射率的格式的结果。在此其中，常用的绝对大气校正方法包括基于辐射传输模型方法（其中包括MOR/DTRAN模型、LOWTRAN模型、ATCOR模型、5S模型、6S模型等）、基于简化辐射传输模型的黑暗像元方法、基于统计学模型的反射率反演方法等；常用的相对大气校正方法包括基于统计的不变目标方法、直方图匹配方法等。针对上述繁多方法的选择，当进行精细的定量遥感研究时，需要采用基于辐射传输模型的大气校正方法；当进行动态监测研究时，可采用相对大气校正或简单的绝对大气校正方法；但所知参数较少时，则只能选择相对简单、对参数要求较少的方法。
  本文采用FLAASH大气校正与QUAC快速大气校正两种方法。
  FLAASH大气校正是一种绝对大气校正方法，其适用于多光谱数据与高光谱卫星影像数据，能够精确补偿大气对辐射的影响。这一方法采用 MODTRAN4+ 辐射传输模型，需要输入地表反射率，通过影像像素光谱中的特征估计大气属性，可以有效去除水蒸气、气溶胶散射等效应带来的干扰，精度较高。而另一方面，通过博客1可知，FLAASH大气校正方法同样对待处理数据有着一定严格的要求，例如需要数据的存储结构为“BIP”或“BIL”模式，像元值类型为经过定标后的辐射亮度（即辐射率）数据，数据单位为(μW)/(cm2nmsr)等。这需要我们针对原始数据加以一定处理。
  QUAC快速大气校正自动从图像中收集不同物质的波谱信息，获取经验值，从而完成高光谱和多光谱的快速大气校正。这一大气校正方法只可以对多光谱、高光谱图像进行处理。
  值得一提的是，FLAASH大气校正与QUAC快速大气校正两种方法得到的结果均为扩大了10000倍的反射率数据。

1.2.4 反演算法

  获得太湖水体区域多个样点的高光谱数据曲线后，分别对其加以一阶微分处理与倒数对数处理。
  本文中，针对太湖水体叶绿素a含量的估计，首先采用常见文献中已给出的模型与参数：

  其中，Chl.a为水体叶绿素a含量，单位为（mg/L）；R_1与R_2为两个特征波段，R_1为近红外波段的反射峰，R_2为红光波段的吸收谷，二者需要通过光谱曲线加以求取；0.325与0.311为两个经验模型参数。一般地，特征光谱曲线其可以依据实测水体叶绿素a含量与光谱数值，通过计算二者相关系数随着对应波长变化而变化的关系，从而加以确定；而经验参数则是依据选择出的特征波段，通过回归分析等拟合手段求出。这一建立模型的正演过程本次我们不做讨论。
  另一方面，由于所获取的高光谱数据可能会受到背景噪声影像，且相似光谱之间因存在明显共线性而导致信息冗余，我们还可以对光谱数据加以一定预处理，从而消除上述噪声，突出光谱特征。往往使用光谱一阶微分、光谱倒数对数与光谱包络线去除等方式加以处理。其中，光谱包络线去除可以方便地在ENVI软件中实现。我们本次采取光谱一阶微分、光谱倒数对数方法。
  光谱一阶微分可以去除部分线性或接近线性的背景，以及噪声光谱对目标光谱的影响。光谱一阶微分计算公式如下：

  其中，R^’ (λ_i )为波段λ_i对应的一阶微分数值，R(λ_i )为波段λ_i原本数值，λ_(i+1)为对应波段中心波长。
  在光谱定量分析中，光谱倒数对数可以十分有效地增强相似光谱之间的差别；计算如下：

  其中，〖R(λ_i )〗_L为波段λ_i对应的倒数对数数值，R(λ_i )为波段λ_i原本数值，λ_i为对应波段中心波长。
  此外，包络线去除法亦是一种常用的光谱分析方法，其可以有效地突出光谱曲线的吸收和反射特征。但我们在此不做讨论，后期我会专门写一篇博客讨论这一内容。
  除水体叶绿素a含量反演外，利用高光谱数据反演地表沉积物颗粒度参数同样被学者们加以广泛研究。有学者借助Hyperion遥感影像，分别建立一个单波段因子与四个多波段组合因子；利用光谱包络线去除方法，挑选出与颗粒度参数相关性较高的几个特征波长，并选用稳定性更高、整体更简单的线性回归模型作为最终的参数反演模型。以平均粒径这一参数为例，具体反演计算模型如下：

  其中，d_m为平均粒径，F为特征波段的运算函数，λ_1与λ_2分别为所选择的两个特征波段。其中，

  同时，本次土壤有机质含量反演模型采用如下公式：

  其中，OM为土壤有机质含量。

2 基于经验比值法、一阶微分法的叶绿素a含量反演

  首先，在未进行大气校正的情况下，我们基于经验比值法、一阶微分法，对叶绿素a含量加以反演。

2.1 数据导入与波段合成

（1） 将数据中“EO1H1190382004095110KZ.tgz”图像压缩包文件解压缩；打开ENVI Classic 5.3（64-bit）软件，选择“File”→“Open Image File”，在弹出的文件选择窗口中分别选中压缩包中后缀名包括8至57的50个“.tif”格式文件；点击“打开”。

（2） 选择→“Layer Stacking”，在弹出的文件选择窗口中选择打开的50个波段图像。

（3） 选择后，需要观察六幅图像在待合成文件列表中的排列顺序。由于后期我们需要对合成后的图像各个波段的信息进行折线图形式的分析操作，因此需要将待合成文件按照“中心波长值”由小至大的顺序排列。点击“Reorder Files”即可实现这一功能。

（4） 结合本文开头提到的三篇博客可知，由于导入“Layer Stacking”模块的波段往往是倒序排列，这使得几十个波段的合成排序操作变得较为复杂——需要进行排序操作百余次。因此本次尝试再将文件导入ENVI软件时就采取倒序导入的方式。但随后发现即使这样操作，也会使得顺序列表中的波段倒序排列。

（5） 顺序排列完毕后，检查投影信息等无误后，点击“OK”即可开始合成操作。

（6） 多次重复这一过程，直至将同一通道内全部波长编号相连的图像各自拼接。

（7） 四个波合成后，将SWIR通道对应的三个一次合成后的图像再一次合成。

2.2 辐射定标与波段合成

（1） 选择“Basic Tools”→“Spectral Math”，在弹出的公式创建窗口中输入本次实验的两个辐射定标公式；输入单个公式完成后，点击下方“Add to List”按钮，即可将键入的公式存入待选择区内。为了减少后期不必要的工作量，可以在编辑完成每一公式后点击“Save”按钮，并将最近一次保存的同一公式文件覆盖，以将待选择区内的两条公式统一保存。若直接选择“Band Math”（如下图）会导致无法对各波段实现简单的批量操作。

（2） 保存公式完成后，点击左下角“OK”按钮，即可开始公式的计算。在弹出的公式变量文件选择窗口中，将这一公式的变量“S1”选择为我们通过上述操作生成、添加的图像文件“layerstacking1_8_57”，并将第二个公式的变量“S1”选择为同样通过上述步骤获得、添加的图像文件“layerstacking5_79_223”。随后，配置对应两个定标输出文件地址等信息。

（3） 配置完成后，点击左下角“OK”按钮，即可开始公式的运行。运行结束，将所得到的研究区两幅辐射定标结果图像导入ENVI软件中，显示如下。

（4） 由两幅图各自的光谱曲线可以看出，V-NIR波段数约为50个，SWIR波段数约为120个。这与我们通过上述波段合成得到的波段数分别为50个、126个的两幅图像相符合。
（5） 再次执行上述波段合成操作，将V-NIR波段图像与SWIR波段图像合成。

（6） 通过本次波段合成，可以看到得到的结果图像包含了V-NIR波段数50个，SWIR波段数126个。

2.3 编辑头文件

（1） 由于在选择波段合成范围时，我们去除了一些“坏波段”，因此合成后的图像中心波长、FWHM等信息均丢失，需要重新将其导入。考虑到176个波段数量较大，我们可以采取Excel数据导出的方式将其填入ENVI影像中。
（2） 首先在Microsoft Office Excel软件中打开原始遥感影像数据文件夹中的“Hyperion_Spectral_Coverage.xls”文件。这一表格文件为我们列出了Hyperion产品Level 1数据242个波段的相关信息，包括波段编号（Hyperion Band）、平均中心波长（Average Wavelength）、FWHM等属性数据。
（3） 将我们最终定标、合成后的176个波段对应的上述三种属性数据导入到一张新的表格中，并导出为TXT文件。其中需要注意，将V-NIR波段与SWIR波段重叠部分的数据剔除。

（4） 其次在“Available Bands List”文件列表中右键选择波段合成后的结果文件，选择“Edit Header”功能。
（5） 选择“Edit Attributes”→“Wavelengths”，将TXT文件导入，并注意将数据单位修改为“Nanometers”。

（6） 此时再次打开上述操作编辑头文件完毕的176个波段的图像，可以看到光谱曲线中横坐标已经由原先的1至176转变为对应波段的中心波长数值范围。

2.4 图像格式转换

经过ENVI处理后，得到的结果图像并不能直接被ERDAS软件所识别，我们需要将其转换为“.img”等格式文件。
（1） 选择“File”→“Save File As”→“ERDAS IMAGINE”，在弹出的文件转换选择窗口中选择经过波段合成与头文件编辑后的结果图像，点击左下角“OK”。

（2） 在弹出的转换文件属性配置窗口中设置，配置好结果图像文件保存路径、保存文件名等。

（3） 文件格式转换完毕后，即可将结果文件在ERDAS软件中打开，并在其中进行后续操作。

2.5 EDRDAS文件导入与裁剪

通过上述步骤得到的“.img”格式文件范围为太湖区域加之其西南部大面积陆地区域，水体占比并不多；而遥感图像整体区域面积较大，文件数据量多；另一方面，我们的主要反演目标区域为太湖水体。因此，我们需要借助自行划定的AOI（Area of Interest）文件，将原本的整体面积地区裁剪为太湖占比较多的地区。
（1） 打开ERDAS IMAGINE 2015软件，在黑色图层窗口区域右键，并选择“Open Raster Layer”，在弹出的文件打开选择窗口中选择经过上述全部预处理后的“.img”结果图像，点击右上角“OK”。ERDAS软件选择文件时，可以借助“Recent”和“Goto”功能，选择最近操作的文件夹路径或文件。
（2） 在图像中划定一个合适的太湖AOI区域，使得全部太湖水体都包含在内。
（3） 选择“Raster”→“Subset & Chip”→“Create Subset Image”，在弹出的文件裁剪配置窗口中选择上述“.img”结果图像，配置输出文件路径和文件名，并选择输出图像文件格式为“Float Single”，在下方“AOI”选项中选择“Viewer”，点击左下角“OK”选项。这样即可以刚刚划定的AOI区域为裁剪范围。

（4） 观察得到的结果图像，若对其直方图加以分析可发现，图像中存在大量数值为零的像元。其实这样并不是图像出现错误，而是在倾斜长方形图像周围区域实际上仍为图像所包括的区域，这些大量黑色的像素数值为0，使得直方图出现零值。

2.6 监督分类

如本文第一部分原理所述，本次实验我们依然利用监督分类方式区分太湖水体、其他水体与其它土地，从而实现更加精确的太湖区域提取。监督分类依然在ERDAS中实现。
（1） 选择“Raster”→“Supervised”→“Supervised Editor”，在弹出的AOI区域显示表中可以看到，此时还没有添加进入任何AOI，表格中处于空白状态。
（2） 分别依据太湖水体、其他水体与其它土地，在图像中划定区域。每一种地表类型需要划出尽可能多的AOI，以期提高最终监督分类的效果与精确度。

（3） 每划定一个AOI，便添加进入“Supervised Editor”中；同一地物类型的AOI添加完毕后，对这一类型的全部AOI加以合并，并赋以一个易于辨认的Value值。

（4） 其中，可以借助卫星地图影像对具体地物类别加以具体区分。但应注意，我们所使用的遥感影像为2004年的数据，而卫星影像数据成像时间相对较为接近我们现在的时间。因此，利用卫星影像亦只可以作为一种参考。
（5） 另一方面，亦发现无论是采用灰度图像还是选择三个波段分别作为R、G、B的值，在“Supervised Editor”中显示的颜色都是一致的，即原有的灰度图像。这一特点与前几次的多光谱博客有所差距，或许为高光谱影像的特征之一。

（6） AOI划分完毕后，保存。选择“Raster”→“Supervised”→“Supervised Classification”，在弹出的监督分类配置窗口中配置输入文件、输出文件与用于指定区域地表类型的AOI文件。

（7） 得到的监督分类结果如下。其中，选择了采取平行六面体规则这一非参数规则。可以看到，由于三种地类分类颜色相近，具体的分类结果需要放大后才可以看清。

（8） 下图中，左图则是采用特征空间规则这一非参数规则加以实现的监督分类结果。二者对比可以看出，右侧的平行六面体规则分类效果明显优于左侧。

2.7 水体光谱曲线提取

（1） 利用“Inquire”模块，分别选取太湖水体中20处不同位置，并上述位置对应的各波段光谱数值导入Excel表格中，作为我们观察特征波段的基础。
（2） 在进行这一步骤时，需要注意在ERDAS软件左下角有一个百分比显示区域。由于176个波段数量较多，数据在复制过程中会有些卡顿；借助百分比显示状态即可看出复制进度。

2.8 特征波段选取与计算

（1） 将20个采样点的光谱曲线导入Excel软件后，利用本文第一部分原理中的公式计算对应光谱数值的一阶微分与倒数对数。并将得到的结果制图。以下三幅图分别为光谱原始数值、一阶微分与倒数对数对应图像。

（2） 依据以上三幅光谱曲线，我们可以从中选取特征波峰或波谷，并将其带入经验公式中。

（3） 如根据上述两幅图，我们分别选择波长为701.5500nm波段与波长为660.8500nm波段。如本文第一部分所示，将其带入文献中所展示的经验比值模型。

（4） 选择“Toolbox”→“Model Maker”→“Model Maker”，在弹出的New_Model窗口中配置模型。
（5） 经过上述ERDAS建模过程，我们得到依据经验比值方法计算的太湖水体叶绿素a含量反演结果图像。可以看到，图像中出现大量负值。这是由于我们直接选用了文献中的参数，且原始图像未经过大气校正造成的。

（6） 此外，依据上述近似的思路，利用光谱的一阶微分数值计算太湖水体叶绿素a含量。其中，本次所选择的一阶微分处理特征波段分别为波长为691.3700nm与671.0200nm的两个波段。

（7） 分别将以上两幅经验比值法、一阶微分法计算得出的叶绿素a含量结果制作为专题地图。上述经验比值法计算得到结果存在较多负值，故此处暂不展示其专题地图——大家继续往后看即可。

3 大气校正及经验比值法波段调整

由以上结果可知，不进行大气校正，所得叶绿素a含量反演结果精度较低，甚至经验比值法计算得到结果存在较多负值，肯定是不对的。因此，这一部分我们基于以下两个方面，对叶绿素a含量反演精度加以提升：

1. 进行大气校正；
2. 对出了问题的经验比值法所选用的波段加以调整。

3.1 转换文件数据格式

（1） 选择“Basic Tools”→“Convert Data (BSQ,BIL,BIP)”，在弹出的文件选择窗口中选择经过波段合成、编辑头文件操作后的结果图像，点击“OK”。

（2） 调整为“BIL”数据存储格式后，即可使用该图像文件进行FLAASH大气校正。

3.2 FLAASH大气校正

（1） 选择“Basic Tools”→“Preprocessing”→“Calibration Utilities”→“FLAASH”，在弹出的转换因子窗口中选择第二项，即单一因子适用于所有波段的情况。由于我们本次所使用数据原有光谱数值为绝对辐射值的标准单位，即(μW)/(cm2nmsr)，这一单位为FLAASH方法所能利用的单位，故我们需要将转换因子设定为1.00。

（2） 随后弹出的配置对话框中，首先选择输入图像文件、输出图像文件目录及名称，同时依据遥感影像的元数据，配置其中心点经纬度、传感器类型（传感器类型一旦选定，系统将会自动匹配传感器高度与像元大小这两个参数）、航行时间；同时依据实际研究区的情况，配置平均海拔高度这一选项；其次，选择合适的地球大气模型和气溶胶模型。其中，地球大气模型需要根据一张标准查找表确定，气溶胶模型则依据研究区域实际情况加以确定。设置FLAASH大气校正参数如下。其中，Hyperion产品Level 1数据的航行时间、图像经纬度等元数据可以通过解压缩后的文件夹中“EO1H1190382004095110KZ_SGS_01.fgdc”文件查阅。

（3） 接下来，选择“Multispectral Settings”。当基本设置中设置了水汽反演模型和气溶胶模型时，我们需要相应地将多光谱相关属性加以配置，配置完毕后点击“OK”。
（4） 全部设置完毕后，可以点击右下角的“Save”按钮，从而将本次对于FLAASH方法的全部属性配置保存；在后期若对同样的遥感数据进行同样的操作时，保存属性配置可以免去多次重复调整相关参数的工作，从而提高效率。
（5） 点击“Apply”，将会开始执行FLAASH大气校正。但在开始前，将会弹出一个FWHM配置窗口。这是由于，FLAASH若针对高光谱图像加以执行，需要配置这一属性。

（6） 通过与导入中心波长相似的操作，我们将FWHM同样通过Excel软件与TXT文件导入。

（7） 随后开始执行FLAASH大气校正。但发现，FLAASH大气校正执行一定时间后，总会出现未包含错误原因的报错提示，并自动暂停执行过程。

（8） 多次尝试依据网络资源对FLAASH大气校正参数加以调整，但这一问题依然存在。针对这一问题的分析附于下文3.4的（7）部分。

3.3 QUAC快速大气校正

（1） FLAASH大气校正无法执行后，尝试使用QUAC快速大气校正方法。选择“Basic Tools”→“Preprocessing”→“Calibration Utilities”→“Quick Atmospheric Correction”，在弹出的文件选择窗口中选择待处理图像文件，点击“OK”。

（2） 得到结果图像文件后，在其上任意位置处右键，选择“Z Profile (Spectrum)”，查看图像中任意位置对应波谱信息。

（3） 将得到的QUAC快速大气校正图像与校正前图像光谱曲线加以对比，可以看到针对地表植被，光谱曲线走势更加符合植被特征曲线，尤其是550nm附近的反射峰、670nm附近的吸收谷、700nm附近的吸收陡坡、900至1300nm附近的平缓趋势等，QUAC大气校正结果图像均有着较好的展现。

3.4 经验比值法调整

（1） 通过ENVI软件QUAC快速大气校正后，尝试将大气校正后的结果图像重新带入第一次未成功的经验比值模型中，再一次计算这种方法得到的叶绿素a含量。

（2） 大气校正后的结果图像为包括大部分非水体的原有范围区域图像，因此需要对其加以重新裁剪、监督分类。

（3） 本次监督分类时，尝试用不同R、G、B波段配置方式。随后发现，若采用Hyperion产品的真彩色配色方式（29:23:16）得到的结果虽然更加符合实际情况，但其图像像素之间变得较为难以识别，不利于监督分类。

（4） 另一方面，在重新进行监督分类时，发现总是会报出如下错误。多次尝试，均无法避免。

（5） 因此，对输入的QUAC大气校正结果图像加以光谱曲线加以复制，并导出到Excel软件中验证。

（6） 可以看到，上述光谱曲线在波长前一部分并无问题，而在930nm左右之后，其数值均不再发生变化。经过多次重复、尝试，找到问题所在——由于波段合成、编辑头文件、大气校正等过程均在ENVI软件中执行，因此使得输出的文件名越来越长（ENVI软件每完成一步光谱处理均会在原波段名称前增添一个操作名称单词）；而过长的ENVI文件名超出了ERDAS软件的文件名识别上限，因此软件将B57波段后的所有波段均识别为同一波段，造成错误。错误的图像如上页最后一幅图；下图则为一幅未经过ENVI大气校正而没有产生这一问题的正常高光谱图像的波段信息。

（7） Excel软件中，可以清晰地看出由B79波段开始，所有光谱数值不再发生变化。之所以在B79波段开始出现这一问题，一定是因其为SWIR波段的第一个波段，而SWIR波段相较之V-NIR波段多经过一次ENVI软件的波段合成（因为其额外需要由三组波段合成为一组）；多经过一次ENVI操作，其名称会多出一个英文单词，从而出现这一问题。且由此分析，FLAASH之所以无法正常完成，亦是因为我是用ENVI进行图像预处理，使得波段名称过长；而FLAASH运行过程中会生成很多新的临时文件，这些新文件的波段命名方式应亦为在各波段前添加单词，从而导致部分波段名称中相互之间唯一不同的字符被“挤出”软件可识别的名称范围，出现重名，导致FLAASH在运行过程中报错。

（8） 因此，为解决经验比值算法出现的叶绿素a含量为负数的问题，决定只得由第二个角度——特征波段入手。因为我们直接运用文献中推导的经验比值模型参数，且特征波段的选择具有一定主观性；因此特征波段亦会影响最终结果取值。

（9） 多次观察反射率光谱曲线，同时适当结合一阶微分、倒数对数光谱曲线 ，以确定出新的特征波段。

（10） 同样利用ERDAS软件，以期建立模型求出采用新特征波段的叶绿素a含量。最终，确定出701.5500nm与691.3700nm这两个特征波段作为选定的波段。上述特征波段确定的叶绿素a含量数值虽仍有小于0的部分，但已有很大改观。

（11） 修正后的经验比值模型求得的太湖水体叶绿素a含量专题地图如图所示（经过指正，以下这幅我的专题地图有误，大家所得结果可能与我的不太一样，依据实际情况判断即可）。同上述第一幅专题地图一致，在制作时需要将图像中太湖周围（即陆地部分）的0值转为NoData。

4 经验比值法、一阶微分法的叶绿素a含量反演结果对比

（1） 针对上述两种水体叶绿素a反演算法得到的两幅专题地图，对比、分析如下。

（2） 其中，左侧图像为经验比值模型反演结果，右侧图像为一阶微分模型反演结果。
（3） 在取值范围上，经验比值结果取值分布于0至0.040mg/L左右，一阶微分结果取值分布于0.020至0.065mg/L左右。依据文献及实际情况，可以看到上述两种方法得到的结果整体均偏小，尤其是缺少含量较高的数值（我国东部亚热带、温带内陆湖泊水体叶绿素a含量最大值往往可达0.110至0.120mg/L以上）；其中，经验比值得到的结果明显小于一阶微分结果，前者最大值甚至仅仅达到后者最小值左右的水平；且经验比值结果中包含极少量负值像素。由此可知，在数值方面，一阶微分模型精度更高一些。当然，我们本次模型中常量参数均直接利用文献给出的数值，而实际中这些参数的数值是需要依据地表水体实测叶绿素含量、各波段与叶绿素含量的相关系数等加以推导、验证得出的，即上述取值范围有很大一部分是由于参数选取并不是很准确导致的。
（4） 在分布趋势上，可以看到两种方法得出的结果均符合“岸边含量达到最高，水体中心含量相对最低，由岸边至水体中心含量逐渐降低”这一情况。而不同的是，经验比值得到的结果随离岸距离的变化极其明显，呈现由东北至西南方向的条带状分布；甚至在太湖这一区域的西南部分，由于河岸向水体中心方向明显延申，导致其随亦是近岸区域，但反演得到的叶绿素a含量并没有随之升高——这样的条带状分布过于明显，可能是由于未经处理的光谱曲线含有的噪声或共线性过多，光谱的大气影响等未很好消除而导致的。而一阶微分得到的结果图像相对前者而言，其叶绿素a含量随水体向内延伸而变化的幅度较小，并不是特别明显。而两幅图中，叶绿素a含量最高和最低的区域则较为一致，均为西南圆弧状岸边处含量整体较高，东北部水域含量整体较低。
（5） 综上所述，可以看出无论是在具体数值取值范围上，还是在数值的空间分布上，一阶微分得到的结果相对较为真实。

5 其他探究内容

5.1 基于ENVI实现微分计算

（1） 在ENVI 5.3 (64-bit) 软件中，可以借助“deriv”函数实现光谱曲线的微分计算。

（2） 与经典版不同的是，ENVI 5.3 (64-bit) 软件在编辑头文件时有所变化，但整体步骤流程一致。此外，太湖区域的提取亦可以在ENVI软件中通过剪裁实现。

5.2 水体表层沉积物平均粒径反演

（1） 依据本文第一部分原理，反演计算太湖区域水体表层沉积物平均粒径。

（2） 求出专题地图。

（3） 由专题地图可以看出，所求得的水体表层沉积物平均粒径取值主要集中于6.0Ф以上；其结果图像中高值分布区域与前期在做监督分类时采用的假彩色图像中太湖水体浑浊（即泛白色）位置（如下图所示）十分一致，可见原本遥感图像中浑浊其实为水体中表层沉积物；沉积物平均粒径越大，水体越浑浊。而上述结果中的条纹状现象可能为特征波段选取并不是很恰当，其间具有部分共线性导致的。

5.3 土壤有机质反演

（1） 依据文献，利用两个特征波段各自倒数对数的一阶微分的比值与土壤有机质含量的二次回归方程，经过同样的建模、制图过程，得出太湖周边地区土壤有机质反演专题地图。
（2） 从中可以看出，土壤有机质含量较高的区域主要集中于城镇中心、道路两旁，这可能是由于人为、交通等导致的碳积累；大部分地区有机质含量处于0.5%至1.0%左右，这个值或许有些小于常见水平。图像中的小面积空白区域多为非太湖湖泊、水田、鱼塘等水体。

本文公式、反演算法等可以参考的文献

欢迎关注公众号：疯狂学习GIS