CN105405112B - 多光谱卫星影像距离偏差指数去雾方法 - Google Patents

Abstract

本发明涉及多光谱卫星影像距离偏差指数去雾方法。根据原始多光谱卫星影像的含雾区域和地标区域的灰度值的差值设定距离偏差指数；然后，将所述距离偏差指数应用在原始光谱卫星影像的局部范围，并对其进行改进。根据上述改进后的距离偏差指数通过计算得出粗糙的雾厚度图像；然后对所述图像进行精细化处理，从原始光谱卫星影像中去除精细化后的雾厚度图像，从而实现对原始多光谱卫星影像的去雾。本发明的方法通用性较强，可以有效减弱原始影像雾区域中的一致性影响，从而恢复雾下方地表的真实光谱信息，获取高质量清晰的地表卫星影像。

Description

多光谱卫星影像距离偏差指数去雾方法

技术领域

本发明涉及多光谱卫星影像距离偏差指数去雾方法，属于光学卫星影像图像处理的技术领域。

背景技术

近地空间环境中，云和雾的影响直接造成光学遥感卫星影像的质量下降和模糊，给影像的大气校正、匹配、定位等预处理和后续的实际应用带来严重阻碍。现有技术中针对卫星图像中云、雾的处理主要通过以下三个方面：(1)基于雾的模糊、变化缓慢等成像特征，利用频率域变换算法消除频率域中的低频成分；(2)利用图像波段间的相关性，在空间域中区分云、雾和真实地表；(3)基于暗像元法，确定原始影像中的低反射率地物，从而分离出雾的影响。

目前现有技术中主要应用的去雾方法为基于暗像通道先验法，这种方法基于计算机视觉中的去雾模型，通过获取透射率图像并估算大气光达到去雾的目的，其主要优势为对雾分布均匀的影像效果显著。但是，当雾霾影响较为严重时，其方法便受到很大限制，主要表现为雾的区域色调突出，出现较大反差。

发明内容

本发明的目的在于克服现有技术的不足，提出了多光谱卫星影像距离偏差指数去雾方法，解决了现有技术中当雾霾影响较为严重或云雾分布不均匀时，去除卫星多光谱影像中的浓雾效果差的技术问题。

本发明是通过如下方案予以实现的：

多光谱卫星影像距离偏差指数去雾方法，步骤如下：

步骤1，根据原始多光谱卫星影像的含雾区域和地标区域的灰度值的差值设定距离偏差指数；然后，将所述距离偏差指数应用在原始光谱卫星影像的局部范围，并对其进行改进；改进后的距离偏差指数表达式为：

其中，在原始光谱卫星图像中，以当前像元点(x,y)为中心，分别在x和y方向上左右上下进行延伸，延伸距离设为line，取该距离内的灰度最小值作为背景值，分别记为min(LineX_i(x,y))和min(LineY_i(x,y))，i表示为波段；DN'_i(x,y)为当雾影响区域作为原始多光谱影像前景像元时，当前像元点(x,y)处的像元灰度值；

步骤2，根据上述改进后的距离偏差指数通过计算得出粗糙的雾厚度图像HTM_i(x,y)；然后对所述图像HTM_i(x,y)进行精细化处理，从原始光谱卫星影像中去除精细化后的雾厚度图像，从而实现对原始多光谱卫星影像的去雾。

进一步的，步骤1中设定距离偏差指数的表达式如下：

其中，为假设正常地表作为影像背景像元时，点(x,y)处的真实灰度值；DN_i'(x,y)为雾影响区域作为影像前景像元时，当前点(x,y)处的像元灰度值；i表示波段；从理论上表示去雾后的像元值，是未知参数。

进一步的，步骤2中利用改进后的距离偏差指数通过计算得出粗糙的雾厚度图像HTM_i(x,y)的表达式如下：

HTM_i(x,y)＝DN_i'(x,y)-DN_i'(x,y)·Δf_i(x,y)

其中，DN_i'(x,y)表示原始影像第i波段点(x,y)处的像素值；DN_i ⁰(x,y)为校正后的像素值；Δf_i(x,y)为距离偏差指数。

进一步的，步骤2中对粗糙的雾厚度图像HTM_i(x,y)的精细化处理具体为：在粗糙的雾厚度图像HTM_i(x,y)中定义大小为ω×ω的不重叠窗口，然后，获取暗像元图像DM_i(x',y')为：

其中，Θ表示以点(x'·ω,y'·ω)为起始点；

然后，对获取的暗像元图像进行中值滤波。

进一步的，对精细化后的雾厚度图像还需要进行修正，进而得到去雾后的原始多光谱卫星影像，则修正去雾后的原始多光谱卫星影像的表达式如下：

其中，minDN_i'为原始多光谱影像中存在不受雾影响的像元。

本发明和现有技术相比的有益效果是：

本发明提出了多光谱卫星影像距离偏差指数去雾方法，本发明将根据含雾区域和地标区域的灰度值定义了距离偏差指数，并对该距离偏差指数进行改进将其应用在多光谱影像的局部范围内，从而获取粗糙的雾厚度图像。然后对该粗糙的雾厚度图像进行精细化处理得到优化的雾厚度图像，从而实现对多光谱影像图像精确的去除。本发明利用定义的距离偏差指数可以有效的从多光谱影像中去除雾厚度图像，尤其当雾霾影响较为严重或者云雾分布不均匀时，其厚度和分布范围变化较大的情况下，可以有效减弱原始影像雾区域中的一致性影响，从而恢复雾下方地表的真实光谱信息，获取高质量清晰的地表卫星影像。并且本发明的方法通用性较强，应用范围较广，对不同类型卫星多光谱影像中的浓雾和阴影具有非常显著的去除效果。

附图说明

图1是本发明实施例的方法流程图；

图2是本发明实施例中原始图像中不同点位的空间分布图；

图3是本发明实施例中原始图像中邻近区域不同点位的距离偏差指数；

图4是本发明实施例中整体和局部的距离偏差指数结果图像；

(a)原始影像；(b)通过公式(7)的结果图像；(c)通过公式(8)的结果图像；

图5是本发明实施例中粗糙的雾厚度图像在经过优化处理得到精细化的雾厚度图像；

(a)粗糙的HTM；(b)优化后的HTM；(c)经滤波后的HTM；

图6是本发明实施例中400×400的Landsat8OLI多光谱影像局部区域的蓝波段雾厚度图像和校正结果；

(a)Landsat8OLI多光谱影像中的原始影像；

(b)Landsat8OLI多光谱影像中精细化后的蓝波段雾厚度图像；

(c)Landsat8OLI多光谱影像中的结果图像；

图7是本发明实施例中400×400的GF-wfv2多光谱影像局部区域的蓝波段雾厚度图像和校正结果；

(a)GF-wfv2多光谱影像中的原始影像；

(b)GF-wfv2多光谱影像中精细化后的蓝波段雾厚度图像；

(c)GF-wfv2多光谱影像中的结果图像；

图8是本发明实施例中400×400的ZY-3多光谱影像局部区域的蓝波段雾厚度图像和校正结果；

(a)ZY-3多光谱影像中的原始影像；

(b)ZY-3多光谱影像中精细化后的蓝波段雾厚度图像；

(c)ZY-3多光谱影像中的结果图像；

图9是本发明实施例的1000×1000的Landsat8OLI影像AHOT和本发明的校正结果及蓝波段直方图；

(a)Landsat8OLI原始影像；

(b)Landsat8OLI影像的AHOT结果；

(c)Landsat8OLI影像在本发明中的结果图像；

(d)Landsat8OLI影像的蓝波段校正前后直方图；

图10是本发明实施例的1000×1000的GF-wfv2影像AHOT和本发明的校正结果及蓝波段直方图；

(a)GF-wfv2原始影像；

(b)GF-wfv2影像的AHOT结果；

(c)GF-wfv2影像在本发明中的结果图像；

(d)GF-wfv2影像的蓝波段校正前后直方图；

图11是本发明实施例的ZY-3影像AHOT和本发明的校正结果及蓝波段直方图；

(a)ZY-3原始影像；

(b)ZY-3影像的AHOT结果；

(c)ZY-3影像在本发明中的结果图像；

(d)ZY-3影像的蓝波段校正前后直方图。

具体实施方式

下面结合附图和实施例对本发明做进一步详细的说明。

多光谱卫星影像距离偏差指数去雾方法，具体步骤如下：

1.对于卫星影像由于其包含复杂的地物，且分布范围广，云和雾的像元灰度值变化在不同地区上是不一致的，因此很难直接通过设定阈值来区分云雾及其下方的正常地表。因此，本实施例中对于原始的多光谱卫星影像，分别从整体图像和局部区域分析原始影像当前像元点(x,y)的距离偏差指数Δf(x,y)变化关系，并计算第i波段的距离偏差指数图像Δf_i(x,y)，通过计算距离偏差指数削弱云、雾的一致性影响因素，从而获取雾厚度图像来区分云、雾及其下方的真实地表；具体方式如下：

如图2所示，在原始的多光谱影像中分布三个点a、b和c，这三个点的灰度值分别表示为DN'_a、DN_b'和DN'_c：设a点处在雾区域中，b和c属于正常区域，由于地表存在高亮度目标，因此雾区域中的a点可能小于正常区域的b和c，所以很难通过灰度值直接对三点进行区分。所以，a、b和c通过如下关系式表示其灰度值：

其中，Ha、Hb和Hc分别为对应点a、b和c的雾的影响，和分别为对应点a、b和c的正常地表的反射能量。

定义，距离偏差指数Δf_i(x,y)为：

其中，为假设正常地表作为影像背景像元时，点(x,y)处的真实灰度值；DN_i'(x,y)为雾影响区域作为影像前景像元时，当前点(x,y)处的像元灰度值；i表示波段。从理论上表示去雾后的像元值，是未知参数。结合式(1)和式(2)，计算点a、b、c三点的距离偏差指数分别为Δf(a)、Δf(b)、Δf(c)，表达式如下：

当点a、b和c三点所受雾的影响相同时，即Ha＝Hb＝Hc，此时的距离偏差指数即可以用来反映三个点真实地表反射能量的比值关系，如下：

如图3所示，在雾的影响区域中，由于雾的变化缓慢，可以近似假设在一定范围内(每一列或每一行)雾的影响是一致的，如图3所示，假设p₁、p₂、p₃三点受雾的影响一致，三者对应的关系为p₁h₁:p₂h₂:p₃h₃，其中p₂和p₃两点较为接近难以区分，而经过距离偏差指数后变为op₁:op₂:op₃，三者间的距离被拉开。因此通过距离偏差指数可以分离这些一致的影响因素，从而复原雾下方真实地表的对比关系。

2.应用距离偏差指数在多光谱影像的整体范围和局部范围

(1)在整体图像范围内应用距离偏差指数：

设原始的多光谱影像的宽高分别为W和H，定义列号和行号分别为j和i，则在列方向上的灰度最小值和最大值分别为X_jMIN和X_jMAX；在行方向上的灰度最小值和最大值分别为Y_iMIN和Y_iMAX，计算影像中单位像元点的灰度值为DN_j,i，其偏离当前列和当前行最小值的距离可表示为：

xdm_j＝DN_j,i-X_jMIN

ydm_i＝DN_j,i-Y_iMIN (5)

根据上述定义的距离偏差指数公式，单位像元点(j,i)在整幅图像上的距离偏差指数可以表示为：

为了使单位像元点在整幅图像中具有相同的地物背景，分别取所有行和列最小值的均值作为固定值代替，即上式可以转化为：

上式(7)即为应用于整体图像的距离偏差指数。

如图4所示，原始影像(图4(a))经过针对整体图像的距离偏差指数运算后得到的结果如(图4(b))所示，雾区域有所收敛，薄雾的影响被减弱，部分地表和阴影区域凸显出来。但是，仍然存在部分缺陷，即图像中存在很多线条，这主要是由于单行单列灰度值最小值不一致所造成的。而且雾厚的地方会出现聚集现象，因此单个像元点相对于整幅图像，其距离偏差指数会随着地表变化呈现区域集中现象(如相邻多行或多列会出现较为一致的块状现象)。基于上述特征，需要多光谱影像局部范围内计算距离偏差指数。

(2)在局部图像范围内应用距离偏差指数：

通过对公式(7)进行改进实现在多光谱影像的局部范围内应用距离偏差指数。即影像的将列和行进行缩小，以影像中的元点(x,y)为中心，分别在x和y方向上左右上下进行延伸，延伸距离设为line，取该距离内的灰度最小值作为背景值，分别记为min(LineX)和min(LineY)，则公式(7)转化为：

公式(8)中得到Δf_i(x,y)即为在多光谱影像的局部范围内应用距离偏差指数。

如图4(c)所示为在局部范围内应用距离偏差指数后得到的结果图像，从中可以看出分布较为均匀的薄雾部分几乎全被消除，阴影和真实地表较在整体范围应用距离偏差指数得到的图像恢复的更为明显，而且保留了图像的整体色调，不足之处是在雾分布较厚的区域及与正常区域相邻的边缘部分，出现雾聚集现象。但是总体而言，根据算法原理应用在局部范围内的距离偏差指数间接反映了正常地表的相互关系，所以本实施例中选取在影像局部范围内应用距离偏差指数的方式用于接下来获取粗糙的雾厚度图像。

3.获取雾厚度图像

由于传感器距离地表较远，所以假设场景对象相对传感器的距离是相同的，所欲传统辐射传输模型如下式所示：

L_λ＝L₀+HR (9)

其中，L_λ表示传感器接收到的辐射亮度，L₀表示地表反射能量和大气程辐射之和，HR则表示雾的贡献值。

又由于原始遥感影像的DN值与辐射亮度具有线性关系，所以上式可化简为式(10)：

其中，DN_i'(x,y)表示原始影像第i波段点(x,y)处的像素值，DN_i ⁰(x,y)为校正后的像素值，HTM_i(x,y)为对应的雾厚度图像。

根据上述定义的距离偏差指数在一定范围内可以反映雾下方真实地表的相互关系，因此，通过局部范围内的距离偏差指数公式(8)可以获取粗糙的雾厚度图像HTM_i(x,y)为：

HTM_i(x,y)＝DN_i'(x,y)-DN_i'(x,y)·Δf_i(x,y) (11)

通过上述公式获取初步的粗糙的雾厚度图像，如图5(a)所示，该图像中除了包含部分地表反射能量，还包括由于算法在x、y方向上取一定距离的单列或单行窗口所导致的十字丝状暗窗口，因此需要对粗糙的雾厚度图像进行进一步的精细化处理。

4.对粗糙的雾厚度图像进行精细化处理

基于上述得到的粗糙的雾厚度图像，为了削减粗糙的雾厚度图像中所包含的地表反射信息，需要对图像进行精细化处理，本实施中采取的精细化处理基本思想为：通过不重叠窗口获取暗像元图像，并对暗像元图像进行滤波处理。具体方法如下：

假设粗糙的雾厚度图像宽高分别为W和H，定义图像中不重叠窗口(即将粗糙的雾厚度图像均分成整数块窗口)大小为ω×ω本实施离中选取ω＝3，则获取的暗像元图像DM_i(x',y')的大小为W/ω和H/ω为：

其中，Θ表示以点(x'·ω,y'·ω)为起始点，图像中不重叠的窗口大小为ω×ω的窗口。

为了消除雾区域和正常区域的边界影响，对暗像元图像采用中值滤波，则滤波后重新得到雾厚度图像如图5(b)所示，针对雾厚度图像部分区域中出现的十字丝状暗窗口，需再进行一次滤波操作，这里同样采用中值滤波，本实施例中窗口大小ω设为9。通过滤波后得到结果如图5(c)。

通过上述滤波后的图像可以在一定程度上削减粗糙的雾厚度图像中包含的地表反射信息。

然后，假定原始影像存在不受雾影响的像元，定义为minDN_i'，因此将细化的雾厚度图像减去对应波段原始影像的灰度值最小值，当差值为负数时不做处理。把最终优化的雾厚度图像定义为HTM_i'(x,y)，进而得到去雾后的原始多光谱卫星影像为：

实验结果与分析：

实验1：

为了验证上述方法的有效性，根据不同影像空间分辨率，分别选取如图Landsat8OLI 30m分辨率河南省登封市区域多光谱影像，时间2013年6月7日；GF-wfv2高分一号宽视场15m分辨率江苏省无锡市区域多光谱影像，时间2013年4月以及ZY-35.8m分辨率北京市某区域多光谱影像，时间2012年11月,以上三幅图像进行实验。对于局部小范围内分布的浓雾或均匀雾，距离偏差指数的延伸距离line设置为15，图像大小均为400×400。图6(a)(b)(c)、图7(a)(b)(c)和图8(a)(b)(c)分别表示Landsat8 OLI、GF-wfv2和ZY-3多光谱影像的原始影像、精细化后的蓝波段雾厚度图像和结果图像。

对于实验图片ZY-3，从图8(a)和8(c)中得出当雾分布比较均匀时，局部区域内雾的贡献值可以认为是相近的，因此距离偏差指数会分离这些一致性影响因素，使图像的整体色调保持一致，达到很好的校正效果；当雾较为严重时，如实验图片Landsat8 OLI和GF-wfv2，此时，距离偏差指数在消除一致性影响因素的同时也会对雾较厚的区域进行削减，还原部分下方地表。此外，对于原始影像中由于对比度产生的阴影，距离偏差指数也具有很好的还原效果，如图6(c)左上角和图7(c)左下角区域。

实验2：

本发明方法不仅可以有效的针对雾霾影响较为严重的情况，而且与其他现有技术的方法相比也具有很显著的效果。在实验2中对于较大范围的区域，由于其所包含的地表种类繁多，此时会有很多高亮度地表，如黄土地等。当延伸距离line较小时，高亮度地表也会出现校正过度现象，因此在处理此类情况时需要对其延伸距离line相应的进行扩大，设置为30，图像大小设置为1000×1000。这里与现有技术AHOT方法进行对比，AHOT方法的基本原理为通过RGB三波段进行雾厚度检测，采用Cloud point校正方法进行雾的去除，下图9到图11分别给出三幅影像的两种方法结果和蓝波段直方图变化。

表1影像质量评价参数表

从图9(b)、图10(b)和图11(b)中可以得出，AHOT算法要求选取不受雾影响的正常区域和较厚的云区域，存在人为的干扰因素，而且要求蓝红波段间具有高度的相关性，其校正结果可以除去大部分雾的影响，但是会使图像的色调发生变化,此外还有一些斑点出现；本发明方法如图9(c)、图10(c)和图11(c)在保留整体色调的同时，对于消除原始影像中的一致性影响因素较为突出。从直方图中可以看出，两种方法的结果图像其蓝波段直方图起始点均与原始影像较相近，即在一定程度上保留图像中的暗像元不变，使得在消除雾的同时，尽可能地保留过多的地表反射能量，防止过度的细节缺失。

表1给出两种方法前后的定量参数指标，其中，表中的original代表原始图像、AHOT代表采用AHOT方法得到的结果图像、Ours代表采用本发明方法得到的结果图像。当雾的影响被减弱时整体图像的均值会有所减小；标准偏差越大则去雾后各像素之间的灰度等级越分散，各像素之间的差异减小；平均梯度越大，图像的清晰度越高，两种校正结果较原始影像均有所提高，本发明方法较AHOT有所偏小，主要原因可能是本文方法在消除雾的一致性影响的同时也削弱了部分地表的反射能量，但是这并不影响地表图像的清晰度。通过多源遥感卫星影像的实验结果证明，本文方法在保存原始图像的色调上较好，具有一定的普适性。

本发明在对糙的雾厚度图像进行精细化处理时通过选取获取暗元图像进行中值滤波的方式，此方法仅为优选的实施例，作为其他方式同样可以采用其他滤波方式的精细化处理，如高斯滤波等。

在本发明给出的思路下，采用对本领域技术人员而言容易想到的方式对上述实施例中的技术手段进行变换、替换、修改，并且起到的作用与本发明中的相应技术手段基本相同、实现的发明目的也基本相同，这样形成的技术方案是对上述实施例进行微调形成的，这种技术方案仍落入本发明的保护范围内。

Claims (3)

1.多光谱卫星影像距离偏差指数去雾方法，其特征在于，步骤如下：

步骤1，根据原始多光谱卫星影像的含雾区域和地表区域的灰度值的差值设定距离偏差指数；设定距离偏差指数的表达式如下：

其中，为假设正常地表作为影像背景像元时，点(x,y)处的真实灰度值；DN_i'(x,y)为雾影响区域作为影像前景像元时，当前点(x,y)处的像元灰度值；i表示波段；从理论上表示去雾后的像元值，是未知参数；然后，将所述距离偏差指数应用在原始多光谱卫星影像的局部范围，并对其进行改进；改进后的距离偏差指数表达式为：

其中，在原始多光谱卫星图像中，以当前像元点(x,y)为中心，分别在x和y方向上左右上下进行延伸，延伸距离设为Line，取该距离内的灰度最小值作为背景值，分别记为min(LineX_i(x,y))和min(LineY_i(x,y))，i表示为波段；DN'_i(x,y)为当雾影响区域作为原始多光谱影像前景像元时，当前像元点(x,y)处的像元灰度值；

步骤2，根据上述改进后的距离偏差指数通过计算得出粗糙的雾厚度图像HTM_i(x,y)；HTM_i(x,y)的表达式为：

HTM_i(x,y)＝DN_i'(x,y)-DN_i'(x,y)·Δf_i(x,y)

其中，DN_i'(x,y)表示原始影像第i波段点(x,y)处的像素值；Δf_i(x,y)为距离偏差指数；

然后对所述图像HTM_i(x,y)进行精细化处理，从原始光谱卫星影像中去除精细化后的雾厚度图像，从而实现对原始多光谱卫星影像的去雾。

2.根据权利要求1所述的多光谱卫星影像距离偏差指数去雾方法，其特征在于，步骤2中对粗糙的雾厚度图像HTM_i(x,y)的精细化处理具体为：在粗糙的雾厚度图像HTM_i(x,y)中定义大小为ω×ω的不重叠窗口，然后，获取暗像元图像DM_i(x',y')为：

其中，Θ表示以点(x'·ω,y'·ω)为起始点；

然后，对获取的暗像元图像进行中值滤波。

3.根据权利要求1所述的多光谱卫星影像距离偏差指数去雾方法，其特征在于，对精细化后的雾厚度图像还需要进行修正，进而得到去雾后的原始多光谱卫星影像，则修正去雾后的原始多光谱卫星影像的表达式如下：

其中，minDN_i'(x,y)为原始多光谱影像中存在不受雾影响的像元。