遥感卫星影像的应用面临一个重要事实：遥感卫星很难获得无云或者雾霾较少的影像，而云或雾霾降低了影像质量，造成影像重要信息的损失，使得依赖于地物表面反射率的科学研究很难开展(Zhang等，2002)。大气污染物、气溶胶、薄云以及地面附近的雾霾等都能影响卫星传感器探测到的光谱辐射值，这在可见光波段表现的尤为显著。对遥感影像进行去雾霾处理是十分必要的，这有助于科学研究和遥感应用，特别是土地覆被的分类和变化检测(Song等，2001)。

目前遥感影像去雾霾的方法大致可以分为两类(Dal Moro和Halounova，2007)。一种是基于辐射传输模型的方法，给定大气参数和光谱辐射值，就可以通过辐射传输模型解算出地物真实反射率。这种方法要求在卫星拍摄影像的时间点前后实地测量的当地大气状况和地表典型地物的表面反射率作为模型的已知条件。这对于很多遥感影像的应用来说难以实施，比如历史演变研究、无人区地表覆被监测等。

另一种是基于单幅影像的雾霾去除方法。最常见的是暗目标减法(Vincent，1972；Chavez，1975)，该方法认为在遥感影像中总是存在一些反射率为零或者很低的像素(称为暗目标)，暗目标的辐射值是由大气作用引起的，将所有像素的辐射值都减去暗目标的辐射值，就可以消除大气的影响。这种方法的局限性在于它只适用于均质条件下的雾霾去除。对于雾霾不均匀(随空间变化)的情况，也产生了许多方法，比如小波变换(Du等，2002)、T-C变换(Richter，1996a，b)和波段组合方法(吴寿江等，2012)等。这些方法大多对地表覆被类型敏感或者只适用于特定的场景，因此限制了它们的适用范围。Zhang等人(2002)提出了一种基于图像变换的方法(Haze optimized transform，HOT)。HOT变换产生的转换图像(HOT图)的振幅与雾霾的影响程度成比例，HOT值的大小反映了雾霾影响程度的强弱。依据HOT值的不同，将原始影像分为多层，每层的大气条件是均质的，然后对每层影像都实施暗目标减法，可以有效消除雾霾的影响。该方法在植被区取得了很好的去雾效果，但是对水体、冰雪覆盖区、裸地表以及人造地物等地表覆被类型却很容易出现过度矫正的情形。Dal Moro和Halounova(2007)针对这种情况提出了改进的方法，其核心思想是为不同类型的敏感地物制作掩膜，然后通过插值计算被去除区域的HOT值。掩膜制作过程中，必须针对不同的敏感地物分别选择阈值，同一地物在不同影像中的阈值也存在差异，这不利于自动化的实现。

本文对Zhang等人(2002)的HOT方法进行了改进。首先，利用原始HOT方法处理Landsat 8 OLI数据，这些数据包含了更多的地表覆被类型和更复杂的组合情况，然后结合实验结果进一步分析HOT的敏感性、过度矫正以及红绿蓝波段RGB合成色彩失真等问题。在此基础上，提出了利用归一化差分植被指数(NDVI)和地物红蓝光谱差(RBSD)制作通用掩膜的改进方法，该方法有效地解决了原始方法存在的地物敏感性缺陷。针对RGB合成的失真，利用百分位数方法确定起始波段的改正值，然后利用散射模型计算其他波段的改正值，统一了不同波段的纠正尺度，避免了某一(或几个)波段的过度矫正，显著改善了RGB合成的效果。

本文使用的Landsat 8 OLI数据下载自USGS网站，空间分辨率为30 m。为便于处理，本文将Landsat 8 OLI的蓝、绿、红、近红外波段合成为.img 格式文件，分别对应于合成影像的1、2、3、4波段。为了充分展示细节，文章中示例影像只是原始影像的部分裁剪。Landsat 8 OLI数据源的详细信息见表 1，表中最后一列表示影像在文章中的引用位置。

利用大气辐射传输模型MODTRAN对Landsat TM数据进行模拟，Zhang等人(2002)发现在晴空和薄雾条件下，Landsat TM数据的1、3(蓝、红)波段高度相关，在2维平面中，可以定义一条晴空线CL(Clear Line)，受雾霾影响时，坐标对(TM1，TM3)在2维平面上的位置会偏移这条晴空线，偏移量与雾霾的影响程度成正比。Zhang等人(2002)定义HOT(如图 1)定量描述偏移量的大小，其表达式如下：

$ {\text{HOT}} = {B_1}\sin \theta- {B_3}\cos \theta $

(1)

式中，B₁、B₃分别代表每个像素1、3波段的传感器记录值(DN值)，θ表示CL与横坐标的夹角。当然，实际的晴空线不一定通过坐标原点，但是平移并不会影响去雾过程。这说明，HOT值的绝对大小并不代表雾霾的绝对厚度；重要的是HOT的相对意义：HOT值越大，雾霾相对越厚。

运用HOT方法去除Landsat TM影像的雾霾和薄云的关键步骤(Zhang等，2002；Zhang和Guindon，2003)包括：

(1) 晴空线定义：手动选择无云区，利用回归模型计算TM3和TM1的相关系数，由此计算出θ值。

(2) 对原始影像执行HOT变换，得到对应的HOT图。

(3) 以HOT图为依据，将原始影像分层，每层影像具有相同的或者某个范围内的HOT值，然后利用直方图计算每层各个波段的最小DN值，进而得到各波段最小DN值与HOT值的关系曲线图(如图 2)。

(4) 以无云区各波段直方图的下界DN值为基准，根据(3)中曲线图，确定不同HOT值对应像素各个波段的改正值，然后实施暗目标减法获得去雾后的影像。例如，假设蓝光波段的基准值为20，对于HOT等于40、蓝光波段最小DN值等于27的单层影像，其蓝光波段的改正值就应该为7。

以Landsat 8 OLI数据为例，上述步骤的实验结果如图 3所示。图 3(a)(b)(c)分别代表植被密集区、农田区和城市区，原图中的红色矩形框为选择的无云区，原图和去雾效果图都是红、绿、蓝3波段合成的真彩色图像。整体而言，HOT图很好地评估了雾霾的强度；去雾后，雾霾的影响被消除，RGB合成影像的对比度得到增强，影像细节变得更加清晰。相比较而言，图 3(a)的去雾效果更好，去雾后的影像颜色与实际相符，未出现过多的偏差。图 3(b)(c)的HOT图多处出现局部极大值和局部极小值，这些区域对应于原图中的水体、裸地、人造地物等地表覆被类型；这说明HOT对这些地物类型敏感，而在去雾后的影像中，这些区域被过度矫正，RGB合成色彩明显失真。

Dal Moro和Halounova(2007)的研究证实，HOT在植被区的应用非常成功，但是诸如水体、冰川积雪、裸地表和人造地物等地表覆被类型则很容易产生过度矫正。他们提出了通过选择阈值为这些类型地物制作掩膜的方法，并将HOT方法发展为两个关键过程：HOT变换和掩膜设计。这样设计的掩膜通用性不够，适用范围受到了限制。

图 3中反映出来的敏感性、过度矫正、RGB合成影像色彩失真等问题极大地限制了HOT的应用，本文将针对这些问题提出相应的改进措施，并通过实验说明改进方法的有效性。

基于前文的分析，可以认为植被覆盖区是HOT图评估雾霾强度的有效区，水体、裸地、人造建筑物等对应的区域则属于HOT的敏感区。因此，可以从原始HOT图中提取有效区对应的像素作为有效像素集，然后根据有效像素集推断敏感区的HOT值；最终得到消除敏感地物影响后的有效HOT图，作为暗目标减法的参考基础。归一化差分植被指数(NDVI)能够有效反映植被信息，因而可以根据V_NDVI确定有效区，并进一步提取有效像素集。NDVI的定义式为：

$ {V_{{\rm{NDVI}}}} =({\rho _{{\rm{NIR}}}} - {\rho _{_{\rm{R}}}})/({\rho _{{\rm{NIR}}}} + {\rho _{\rm{R}}}) $

(2)

式中，ρ表示地表反射率，NIR和R分别代表近红外波段和红波段。严格来说，这里的反射率为消除大气影响后的真实反射率，考虑到在雾霾去除之前无法根据原始影像获取真实反射率，将直接利用像素的DN值进行计算。因此，本文的NDVI计算式为：

$ V_{\text{NDVI}} =(DN_{\text{NIR}} - DN_{\rm{R}})/(DN_{\text{NIR}} + DN_{\rm{R}}) $

(3)

V_NDVI的取值范围为[-1，1]，负值表示地面覆盖为云、水、雪等；0表示有岩石或裸土等；正值表示有植被覆盖，正好对应于HOT图评估雾霾强度的有效区。由于计算V_NDVI过程中未排除大气作用的影响，本文中提取植被过程中选择的阈值将会小于0(经验表明，-0.1较为合适)。同时需要指出的是，厚云覆盖区域因为地物信息被遮盖，即使存在植被也会表现出较低的V_NDVI值，这些区域也一并被排除在有效区外。即便对这些区域去雾，恢复出来的信息也是不真实的。

通过实验发现，地表偏蓝的非植被地物类型也表现出较高的NDVI指数，这对于根据NDVI提取有效区非常不利。根据蓝色地物蓝光波段反射率高、红光波段值反射率低的特点，可以计算地物红蓝光谱差值RBSD(Red-Blue Spectral Difference)，这里依然使用像素DN值代替地表反射率，其计算式为：

$ V_{\text{RBSD}} = DN_{\text{blue}} - DN_{\text{red}} $

(4)

V_RBSD异常大的像素就是需要检测的蓝色地物，对于异常大值的临界点，可以通过V_RBSD的直方图确定。此外，对于地表呈现红色、棕色(一般为裸地表)的地物，它们的V_RBSD值表现为负值，而植被区则为相对较低的正值，这一特点可作为NDVI提取植被区的补充。图 4展示了利用V_RBSD检测高蓝地物的效果，图 4(b)中白色表示蓝色地物，黑色表示红色或者棕色地物，图中左下角的白色为原图中的水系，可见RBSD同时实现了敏感水体的检测。

以NDVI和RBSD为依据，确定有效区；用不同的标签分别标记有效区和敏感区，就可以设计出理想的通用掩膜。NDVI区分植被和非植被的阈值相对固定，利于自动化的实现；并且可以一次性排除多种敏感地物，简化了制作掩膜的过程，制作的掩膜通用性更强。假设V_NDVI的阈值为T₁，V_RBSD的阈值为T₂₁(下临界值)、T₂₂(上临界值)，H表示大小等于原始影像的掩膜矩阵，H(i，j)表示行号为i列号为j的像素，用0标记敏感区，1标记有效区，则掩膜可描述如下：

如果V_NDVI＞T₁并且T₂₁＜V_RBSD＜T₂₂，则H(i，j)=1；否则，H(i，j)=0;

将H作用于原始HOT图得到有效像素集，然后通过插值从有效像素集推算敏感区的HOT值，最终生成有效HOT图。本文利用3×3模板处理有效像素集，将模板内有效像素点的HOT平均值作为模板中心像素点的HOT值。

Zhang等人(2002)的分析证明，在雾霾和薄云环境下，HOT对大气的物质组成类型不敏感，利用掩膜去除敏感地物后，地表类型的影响也被消除，因此可以认为，HOT值的大小只与雾霾强度相关，即HOT值有效地反映了雾霾的影响程度。参照Zhang等人(2002)的思路，本文以有效HOT图为参考，对欲去雾的波段实施分层，每层影像具有相同的或者某个范围内的HOT值(这由分层的间隔决定)，分层的间隔可以自由选择，一般选择间隔值为1—5。分层过后，每层影像的大气条件都是均质的，满足了暗目标减法的前提条件。因此，对每层影像分别实施暗目标减法，就能有效消除雾霾的影响。

Zhang等人(2002)在实施暗目标减法的过程中，相互独立地寻找每个波段的暗目标；本文把同一层不同波段的暗目标称为一个暗目标组。这样确定的暗目标组一般来自不同的像元，这些像元的大气状况存在差异，致使暗目标组中不同波段的值不符合真实的大气散射模型。从而导致影像某一(或几个)波段被过度矫正，去雾后影像的波段间关系不符合相对大气散射模型，出现RGB合成影像颜色失真问题。为解决此问题，本文并不独立地寻找各个波段的暗目标，而是考虑大气散射模型；去雾过程中，选择某一波段(一般为蓝光波段)作为起始波段，寻找该波段各个层对应的暗目标，并确定该波段各层的改正值大小，而其余波段的改正值将由波段间散射强度的相互关系直接产生。

一般而言，散射光的强度与入射光波长成反比，Chavez(1988)认为在雾霾环境下这种反比关系可以限制在λ^-0.7与λ^-0.5之间。本文选择λ^-0.7关系进行实验，表 2列出了根据该关系计算的Landsat 8 OLI数据2—7波段的散射强度值。表 3列出了分别以Landsat 8 OLI 2、3、4(蓝、绿、红)波段为起始波段计算其余波段改正值时的乘算因子。以蓝光波段作为起始波段为例，如果蓝光波段的改正值为AD_blue，则绿、红、近红外波段的改正值分别为：

$A{D_{{\text{green}}}} = 0.898 \times AD_{{\text{blue}}}$

$A{D_{{\text{red}}}} = 0.807 \times AD_{{\text{blue}}}；$

$A{D_{{\text{nir}}}} = 0.664 \times AD_{{\text{blue}}}。$

对起始波段(蓝光)进行分层处理后，按本文第3部分步骤(3)、(4)处理，可以获取上述AD_blue值。但是利用直方图最小下界确定的暗目标值通常不稳定，这在蓝光波段表现的尤为显著(如图 2蓝色线)。在去雾过程中，相邻层间的改正值跳跃性改变，去雾后的影像会出现许多块斑和光晕(如图 5(a))。解决此问题的一种有效方式是选择各个直方图的百分位数作为替代，由图 2可以看到25/50/75th 百分位数线(Percentile)表现得比(Lower Bound)更加稳定，由此确定的相邻层间的改正值将平滑地改变，去雾后的影像局部整体性更强，不存在块斑和光晕(如图 5(b))。

采用百分位数方法后，基准值的确定方法也相应地发生改变。从图 2可以看出，直方图的百分位数曲线先单调递减后单调递增，Dal Moro和Halounova(2007)的实验也体现出相似的规律。百分位数曲线代表不同层暗目标的DN值，那么该曲线的最小值则对应于全幅影像的暗目标，这正是暗目标减法过程的基准值。本文将此最小DN值作为后续暗目标减法过程的基准。同时，最小DN值在曲线中对应的HOT值也确定了雾霾或薄云的HOT范围的下界值。

在Zhang等人(2002)的方法中，当影像存在厚云的时候，暗目标减法过程中还需要确定雾霾或薄云HOT范围的上界值。本文4.1节利用NDVI设计通用掩膜的过程中将厚云部分划为敏感区，经过4.3节的掩膜处理后，在有效HOT图中厚云对应的HOT值由其最近邻的薄云或无云区的HOT值插值得到，插值后的值远远小于原始的HOT值。本文在暗目标减法过程采用相同的方式处理厚云和薄云，薄云的影响被消除，而厚云部分得到保留，区别在于去雾后厚云区的DN值比原始值要略微减小，但这并不妨碍厚云的识别和提取(如图 6)。

改进后的HOT方法，具体操作流程如图 7所示。原始HOT方法对地表覆被类型敏感，大部分地区容易产生过度矫正，去雾后的RGB色彩合成效果不佳。鉴于HOT方法在植被覆盖密集区较好的去雾效果，改进的HOT方法将植被覆盖区作为HOT评估雾霾强度的有效区，有效区对应的HOT图部分为有效像素集。然后根据有效像素集去推断敏感区的HOT值，得到有效HOT图。最后以有效HOT图为参考，实施暗目标减法，消除雾霾的影响。为了提取有效区，需要制作掩膜。归一化差分植被指数(NDVI)可以很好地检测植被，但是对高蓝地物敏感，而地物红蓝光谱差(RBSD)能检测高蓝地物，二者相结合实现了通用掩膜的制作，并且避免了水体、裸地、人造地物、冰雪等敏感地物的干扰。在暗目标减法过程中，本文利用取直方图百分位数的方法确定起始波段的改正值，利用散射模型计算其余波段的改正值，统一了不同波段间的纠正尺度，解决了某一(或几个)波段过度矫正的问题，避免了块斑和光晕的产生。同时，改进后的方法可以区分影像中的厚云和薄云，并消除薄云的影响。

图 3最右一列展示了改进后的效果图；图 8列出了实验原图以及改进前后的HOT方法去雾效果图；图 9展示了图 8(a)、图 8(b)与图 8(c)的局部效果对比，它们的具体位置对应于图 8(a)中带有数字标记的矩形框，涵盖的土地覆被类型包括：植被(a)、水体(b)、城市建筑区(c、d)以及裸地表(e)；图 10展示了地表覆被存在冰雪的情况。对比分析可以看出，改进后的HOT方法能够有效消除雾霾和薄云的影响、提高影像清晰度；由于DN值减小，改进后的去雾影像亮度降低。对比改进前后的去雾效果，改进后的方法有效避免了水体、裸地表、人造地物以及冰雪覆盖区等产生的敏感性问题；各类典型地物的色彩和结构与实际情况相符，有效解决了原始HOT方法产生的过度矫正、RGB合成色彩失真问题。

图 11的散点图更直观地反应了去雾后红蓝光谱空间的变化。图中第一行的样本取自图 8(a)中矩形框C和矩形框1，第二行的样本取自图 8(a)中矩形框C和矩形框3，C代表无云区，1代表植被覆盖密集的雾霾区，3代表地人造地物密集的雾霾区。图中左列中轴线标红的部分反应了无云区的光谱空间，其轴线即Zhang等人(2002)定义的晴空线(CL)；区域1、3受到雾霾的影响，它们的中轴线相对晴空线产生偏移，去雾后中轴线与晴空线基本吻合。对比可以发现，第二行右列相比中列保留了中轴线右下侧部分(图中红色矩形框标记部分)，这说明改进后的方法有效避免区域3中人造地物的过度矫正(如图 9(c))。

对于图 9(f)，去雾后的影像明显偏白，类似的情况还出现在植被覆盖区内的裸地、小区域人工建筑区，如图 8(c)左上角区域。这些区域在暗目标减法过程中，减掉的值比正常值偏低，这主要是由插值过程造成的。这些区域内部本身缺少植被，在执行掩膜操作后，其自身的HOT信息消失，插值后的HOT值由周边植被推算而来；另一种可能是该区域保留的HOT值对于插值点的影响弱于周边植被。由图 8(a)可见，周边的植被为相对的无云区，从而致使插值后获得的HOT值比原本的值要低，在去雾过程中确定的改正量相应地也就比正常值低。这是本文改进HOT方法的一个缺陷，该方法必须依赖于区域内的植被，植被缺失将导致该区域的去雾结果产生偏差。

总结起来，本文的改进包括以下几个方面：

(1) NDVI结合RBSD制作通用的HOT掩膜

NDVI区分植被和非植被的阈值相对固定，利于自动化的实现；并且可以一次性排除多种敏感地物，简化了制作掩膜的过程，制作的掩膜通用性更强。

(2) 使用百分位数确定起始波段改正值

直方图的百分位数相比直方图下界更加稳定，能够保证不同HOT层间的改正值相对平滑地变化，避免块斑和光晕的出现。

(3) 考虑散射模型

利用散射模型从起始波段推算其余波段的改正量的大小，保证不同波段的改正量符合相对大气散射模型，使得RGB合成后的色彩与实际相符。同时也减少了逐波段寻找暗目标的时间消耗。

当然，改进后的方法也还存在许多方面需要进一步完善：

(1) 该方法要求以植被作为基础，对于植被分布少的区域(如西北干旱区)、植被未生长季节(如北方冬季)的影像很难处理。同时，对于图 9(f)反映出来的情况难以有效解决。

(2) 掩膜制作方法依然会漏掉某些特殊的敏感地物，而且掩膜过程去除的部分也并不完全属于异常值。另一个思路是直接检测异常值，然后用正常区域的HOT值去填充异常区域。但是，异常值的自动检测和填充是这一思路的难点，也是未来研究的重点。

(3) 缺乏客观的评价体系。如何建立评价指标或者评价函数来衡量去雾后影像的质量也是未来研究的一个方向。