CN104966273A - 适用于光学遥感影像的dcm-htm去雾霾方法 - Google Patents

Abstract

本发明具体涉及一种适用于光学遥感影像的DCM-HTM去雾霾方法，包括以下步骤，步骤一，雾区检测和提取；步骤二，雾厚度图像HTM提取；对获取的暗通道图像选取大小为ω^r×ω^t的不重叠窗口HTM_W(ω^t，ω^t)，ω^t∈[3，9]求取雾厚度图像HTM；步骤三，雾区校正；通过雾厚度图像HTM和所提取的fog_region二值图像，求取地物背景贡献值BACK，以此来修正HTM，最后基于辐射传输原理，在检测出的DCM雾区内，通过HTM去雾模型对原始影像进行校正处理。本发明通过对雾区域范围内的非亮度目标取小窗口实现在保证高亮度目标不被过度校正的基础上达到进一步削减雾的目的；基于检测的雾区域，结合雾厚度图像HTM计算正常区域地物背景贡献值BACK来修正HTM，很好的保留了正常区域的灰度值。

Description

适用于光学遥感影像的DCM-HTM去雾霾方法

技术领域

本发明涉及光学遥感影像领域，尤其涉及一种光学遥感影像的去雾霾方法。

背景技术

受到云、雾霾等大气环境因素的影响，光学卫星遥感影像的成像质量下降，影像的判读和量测性能降低，制约了影像的进一步应用。如以2014年度USGS官网发布的包含登封检校场区域的同一视场多时相Landsat8 OLI 30米分辨率光学卫星影像为例，全年受云影响较为严重的有4个月份，分别是4、6、7和12月；而受雾霾影响的主要有五个月份，分别是4、5、6、7和11月，占了全年的50％。因此，如何有效消除或减小云和雾霾对光学遥感影像的影响，充分发挥影像的效能，已成为对地观测领域所关注的一个热点问题。

从遥感物理上来讲，薄云和雾霾均分布在近地空间，对可见光波段影像的散射作用基本一致的，在卫星影像上的成像特征较为相似，常表现为具有较高亮度、模糊等特征。所以，大部分去雾方法将二者归为一类进行处理，为表述方便，本文也将薄云和雾霾统称为雾。

在计算机视觉领域，比较典型的去雾方法有暗通道先验理论(Srinvasa and Shree，2002，2003；何凯明等，2011)、同态滤波(赵忠明和朱重光，1996)和Retinex(Edwin，1986)以及在以上基础上进行的改进等(庞嘉豪等，2011；何凯明等，2010；兰霞等，2013；何仁杰等，2013；郭璠等，2014)，针对的主要是近距离摄取的风景图像，特点是在成像范围内雾是均匀分布的。因此，很多学者也将此类方法应用于雾分布均匀的遥感图像校正，取得了不错的效果(龙焦等，2014；吴小平等，2012；胡长苗，2014)；但是，一景光学卫星遥感影像的地面覆盖范围往往比较大，包含的地面特征繁杂，雾的分布通常是不均匀的。所以，当雾厚度不一、分布不均时，现有的去雾方法大多收效甚微，且处理后不可避免的会影响原来正常区域的成像质量。

针对此类问题，很多学者在受影响区域的检测和校正方面做了很多研究。如在云覆盖检测评估方面，Landsat长期收购计划LTAP(Long Term Acquisition Plan)通过采集全球多个样点数据对Landsat7和ETM+，推出了自动云检测评估系统ACCA(the Automated Cloud-Cover Assessment)，对薄云、厚云和正常区域进行分类(Hollingsworth等，1996)；由于该算法多次用到了亮度温度图像，(Pasquale等，2012)针对Landstat8 OLI(the Operational Land Imager)数据，通过合成亮度温度图像，发展了AT-ACCA(Artificial Thermal-Automated Cloud-Cover Assessment)算法，该算法尽管可以很精确地估算出云的覆盖区域，但也只限于Landsat系列影像。(Aliaksei等，2014)通过改进暗目标提取的方法计算雾厚度图像(HTM)，并分别对Landsat 8 OLI和WorldView-2进行雾区域检测和校正，去雾后的图像其光谱特征变化具有很好的一致性。算法主要的不足之处在于处理沙漠、雪地等高反射背景图像时，难以通过获取暗目标来估算HTM，此外，其方法需要进行波段合成，并不适用于缺乏短波段的影像。

如上所述，现有的去雾图像处理算法或只针对雾全局均匀分布的图像，或只能处理某一类型影像，或难以对某些地物特征的影像进行有效处理，算法的通用性较差。

发明内容

本发明的目的在于提供一种结合暗通道图像和雾厚度图像的DCM-HTM算法，在对光学卫星遥感影像雾区进行精确检测的基础上，实现了可针对局部雾区的有效去雾校正。

为实现上述目的，本发明提供一种适用于光学遥感影像的DCM-HTM去雾霾方法，包括以下几个步骤：

步骤一，雾区检测和提取；通过选取大小为ω×ω的暗通道窗口darkC_W(ω×ω)，对宽高为(W，H)的原始影像求取暗通道图像DCM，然后对DCM进行mean-shift平滑并根据直方图自动设定阈值生成二值图像fog_region，二值图像fog_region中灰度值为1的区域即为雾区域；

步骤二，雾厚度图像HTM提取；对获取的暗通道图像选取大小为ω′×ω′的不重叠窗口HTM_W(ω′，ω′)，ω′∈[3，9]求取雾厚度图像HTM，由于通过不重叠窗口后所获取的HTM是原始影像宽高的1/ω′，对其再进行滤波、重采样回原图大小；

步骤三，雾区校正；通过雾厚度图像HTM和所提取的fog_region二值图像，求取地物背景贡献值BACK，以此来修正HTM，达到保留正常区域的目的，最后基于辐射传输原理，在检测出的DCM雾区内，通过HTM去雾模型对原始影像进行校正处理。

其中，在所述步骤一中，暗通道图像I^dark(x)的获取方法如公式 所示；其中，I^c表示原始图像的每个通道，c表示R、G、B三通道，Ω(x)表示以像素x为中心的一个窗口；

定义I_i表示第i波段原始影像，DCM_i(x)为其对应的暗通道图像，其获取方法如公式所示；

定义总的概率密度为两个规则概率密度和乘积的结果，其中，x³是坐标空间部分，x^r是对应的颜色空间部分，h_s和h_r分别代表坐标空间和颜色空间的核函数带宽，c为归一化的相关系数，用于滤波并分割时所需设定的变量用h＝h(h_s，h_r，M)表示，其结果将包含像素点少于数值M的区域与它最相似的区域合并。

其中，在所述步骤一中，经过mean-shift平滑处理后，利用公式  ${\mathrm{DN}}_{\mathrm{mean}\_\mathrm{shift}}^{\&prime;}=\left\{\begin{array}{cc}0,& {\mathrm{DN}}_{\mathrm{mean}\_\mathrm{shift}}<T\\ 1,& {\mathrm{DN}}_{\mathrm{mean}\_\mathrm{shift}}>T\end{array}\right.,$ 根据阈值T对暗通道图像二值化，将图像分为背景区域(DN′_{mean_shift}＝0)和雾区域(DN′_{mean_shift}＝1)；在^[0，100]区间内直方图最大值i_histMax和最小值i_histMin处的DN值DN[i_histMax]和DN[i_histMin]来确定最佳阈值T，具体求值方法如公式T＝(DN[i_histMax]+DN[i_histMin])/2DN[i_histMax]，DN[i_histMin]≠0所示。

其中，在所述步骤二中，首先假定暗通道图像DCM宽高分别为W和H，根据所选取的不重叠窗口HTM_W(ω′，ω′)计算初始输出的HTM大小，即为DCM宽高的1/ω′；则从输出的HTM像元坐标点(x′，y′)出发，读取DCM起始点(ω′*x′，ω′*y′)，宽高均为ω′大小的窗口像元值，记为Ω(x′，y′)，最后取该窗口内的最小值作为当前点的HTM值，计算方法如公式HTM_i(x′，y′)＝min_{y∈Ω(x′，y′)}{DCM_i(y)}所示；其中，DCM_i(y)为第i波段暗通道图像在窗口Ω(x′，y′)范围内的像元值，HTM_i(x′，y′)则表示对应的宽高为(W/ω′，H/ω′)的雾厚度图像。

其中，对HTM进行滤波，选择中值滤波，模板大小为m×m，m∈[3，9]，来平滑矩阵窗口所产生的噪声并重采样成原图大小。

其中，所述HTM包括三部分，分别是大气中气溶胶等的散射部分、一些高亮度地表以及雾的影响部分；这里将前两者归为大气的影响，以背景BACK表示，则修正后的雾厚度图像HTM′用 ${\mathrm{HTM}}^{\&prime;}=\left\{\begin{array}{cc}0& ,\mathrm{HTM}\&le;\mathrm{BACK}\\ \mathrm{HTM}-\mathrm{BACK}& ,\mathrm{HTM}>\mathrm{BACK}\end{array}\right.$ 表示；则 其中，DN′_i(x，y)为原始影像第i波段点(x，y)处的灰度值，为其校正的结果图像，HTM_i(x，y)为雾厚度图像。

本发明的有益效果是：本发明提供的适用于光学遥感影像的DCM-HTM去雾霾方法，基于暗通道图像包含少量的地表反射信息这一优势条件，采用mean-shift和直方图自动阈值实现较为精确地雾区域检测和提取；通过对雾区域范围内的非亮度目标提取小窗口实现在保证高亮度目标不被过度校正的基础上达到进一步削减雾的目的；基于检测的雾区域，结合雾厚度图像HTM计算正常区域地物背景贡献值BACK来修正HTM，很好的保留了正常区域的灰度值。

附图说明

图1为本发明提供的适用于光学遥感影像的DCM-HTM去雾霾方法的流程图；

图2为实验1中的原始影像图片数据1；

图3为实验1中对数据1进行DCM-HTM校正后的结果图片；

图4为实验1中的数据1校正前后灰度值变化图；

图5为实验1中的原始影像图片数据2；

图6为实验1中对数据2进行DCM-HTM校正后的结果图片；

图7为实验1中的数据2校正前后灰度值变化图；

图8为Landsat80LI的蓝波段的高亮度目标图像；

图9为Landsat80LI的实验1暗通道图像；

图10为Landsat80LI的实验2暗通道图像；

图11为Landsat80LI未考虑ki校正结果的(R、G、B)显示；

图12为Landsat80LI考虑ki校正结果的(R、G、B)显示；

图13为GF-1wfv2的蓝波段的高亮度目标图像

图14为GF-1wfv2的实验1暗通道图像；

图15为GF-1wfv2的实验2暗通道图像；

图16为GF-1wfv2未考虑ki校正结果的(R、G、B)显示；

图17为GF-1wfv2考虑ki校正结果的(R、G、B)显示；

图18为Landsat80LI校正前后图像的直方图变化图；

图19为GF-1wfv2校正前后图像的直方图变化图。

具体实施方式

参阅图1，本发明提供的适用于光学遥感影像的DCM-HTM去雾霾方法，包括以下几个步骤：

步骤一，雾区检测和提取；通过选取大小为ω×ω的暗通道窗口darkC_W(ω×ω)，对宽高为(W，H)的原始影像求取暗通道图像DCM，然后对DCM进行mean-shift平滑并根据直方图自动设定阈值生成二值图像fog_region，二值图像fog_region中灰度值为1的区域即为雾区域；

步骤二，雾厚度图像HTM提取；对获取的暗通道图像选取大小为ω′×ω′的不重叠窗口HTM_W(ω′，ω′)，ω′∈[3，9]求取雾厚度图像HTM，由于通过不重叠窗口后所获取的HTM是原始影像宽高的1/ω′，对其再进行滤波、重采样回原图大小；

步骤三，雾区校正；通过雾厚度图像HTM和所提取的fog_region二值图像，求取地物背景贡献值BACK，以此来修正HTM，达到保留正常区域的目的，最后基于辐射传输原理，在检测出的DCM雾区内，通过HTM去雾模型对原始影像进行校正处理。

相较于现有技术，本发明提供的适用于光学遥感影像的DCM-HTM去雾霾方法，基于暗通道图像包含少量的地表反射信息这一优势条件，采用mean-shift和直方图自动阈值实现较为精确地雾区域检测和提取；通过对雾区域范围内的非亮度目标取小窗口实现在保证高亮度目标不被过度校正的基础上达到进一步削减雾的目的；基于检测的雾区域，结合雾厚度图像HTM计算正常区域地物背景贡献值BACK来修正HTM，很好的保留了正常区域的灰度值。

结合两个实验，具体说明。

实验一，选取了两组数据开展实验，数据1是河南省登封区域Landsat80LI30m分辨率多光谱卫星影像数据，影像获取时间2014年6月7日，影像大小为(1000×1000)，数据2是江苏省无锡市鄱阳湖区域16m分辨率的宽视场GF-1wfv2影像，获取时间2013年8月9日，影像大小为570×690。暗通道窗口均取(15×15)，不重叠窗口大小取(3×3)。

实验从原始影像出发，依次获取暗通道图像DCM、雾厚度图像HTM和校正结果图像。Landsat80LI的校正处理结果如图2-4所示，图中给出原始影像和校正结果并以(R、G、B)彩色图像进行显示。

高分宽视场GF-1wfv2影像数据校正结果如图5-7所示。同样给出原始影像和校正结果并以(R、G、B)彩色图像进行显示。

分别从图4和图7的校正前后图像中任取相对应的一行数据，当HTM低于背景BACK值时，不做处理，相当于保留正常区域的原始灰度值。

从图3的校正结果来看，在原图左上方的山区和右下方城市区域，雾校 正模型对均匀薄雾的效果比较明显，很好的恢复其下方的真实地表；但是在浓度分布不均的地方仍然存在厚度不等的残留，其主要原因是雾厚度图像HTM的理想来源应是地表反射作用较小而雾影响较大的短波段原始影像，因此随着波段的增长，通过不重叠窗口所获取的HTM越不够精确；从图6可看出其原始影像受雾霾影响非常严重，校正过后，改善效果也极其明显。

由于mean-shift可以使图像中的高亮度地物凸显出来，而分割结果有时会将诸如旱田等中高亮度反射的地物遗漏。基于上述现象，通过对原始图像进行mean-shift均值平滑，并设定阈值(取平滑图像的均值)来确定高亮度目标bright_target(如高反射房顶、旱田等)，所以定义k_i为第i波段在给定原始暗通道窗口大小的前提下，高亮度地表区域的窗口变化条件，如公式

$k_i=\left\{\begin{array}{cc}\mathrm{wndSize}/2,& \mathrm{bright}\_t\arg \mathrm{et}\&NotElement;\mathrm{fog}\_\mathrm{region}\\ \mathrm{wndSize},& \mathrm{else}\end{array}\right..$

其中，wndSize为初始暗通道窗口大小，fog_region和bright_target分别为检测的雾区域和高亮度目标，二者均已二值化，即fog_region＞0，bright_target＞0时分别表示为雾区域和高亮度目标；所以当高亮度目标在雾区中时，为了不使高亮度目标出现校正过度现象，窗口大小不变；当高亮度目标不在雾区时，雾区窗口设为原来的一半，以尽可能多的削减雾的影响，而正常区域则以原窗口大小。

实验二

实验二基于公式 $k_i=\left\{\begin{array}{cc}\mathrm{wndSize}/2,& \mathrm{bright}\_t\arg \mathrm{et}\&NotElement;\mathrm{fog}\_\mathrm{region}\\ \mathrm{wndSize},& \mathrm{else}\end{array}\right.$ 给出的条件分别对前面实验1数据进行校正，并与实验1校正结果进行比对，暗通道窗口设置均设为(15×15)，不重叠窗口设为(3×3)。

Landsat80LI的校正处理结果如图8-12所示，这里给出蓝波段的高亮度目标图像、实验一暗通道图像和对应的实验二暗通道图像，并与实验一局部校正结果进行(R、G、B)显示。

高分宽视场GF-1wfv2影像数据校正结果如图13-17所示，同样给出蓝波段的高亮度目标图像、实验一暗通道图像和对应的实验二暗通道图像，并与实验一局部校正结果进行(R、G、B)显示。

从图7-17所示的校正结果可以看出，原始图像中均分布着厚度不均的雾和霾，给原始影像造成了严重的降质影响，通过考虑k_i条件自动选择暗通道窗口，图像保真度基本不变，相比于未考虑时所消除的雾要多一些，因此地表的真实颜色得以凸显。

图18和图19给出了实验一实验二两组实验的直方图，可以看出，两组实验的直方图较原始影像均有明显的偏移，而且起始增长曲线与原始影像较为一致，间接说明了算法在保留正常区域上有明显的效果。

表1给出两组实验前后相关的质量参数指标。从中可以看出，两组实验结果较为接近，均值部分变化反映了原始影像与雾厚度图像间的差值关系，雾厚度图像值越大，原始雾区域削减的越多；其峰值性噪比和平均梯度较原始数据均有明显提高，表明原始影像经过校正后图像的清晰度提高，可辨识性增强。

下表给出两组实验前后相关的质量参数指标。从中可以看出，两组实验结果较为接近，均值部分变化反映了原始影像与雾厚度图像间的差值关系，雾厚度图像值越大，原始雾区域削减的越多；其峰值性噪比和平均梯度较原始数据均有明显提高，表明原始影像经过校正后图像的清晰度提高，可辨识性增强。

以上仅为本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种适用于光学遥感影像的DCM-HTM去雾霾方法，其特征在于，包括以下几个步骤：

步骤一，雾区检测和提取；通过选取大小为ω×ω的暗通道窗口darkC_W(ω×ω)，对宽高为(W，H)的原始影像求取暗通道图像DCM，然后对DCM进行mean-shift平滑并根据直方图自动设定阈值生成二值图像fog_region，二值图像fog_region中灰度值为1的区域即为雾区域；

步骤二，雾厚度图像HTM提取；对获取的暗通道图像选取大小为ω′×ω′的不重叠窗口HTM_W(ω′，ω′)，ω′∈[3，9]求取雾厚度图像HTM，由于通过不重叠窗口后所获取的HTM是原始影像宽高的1/ω′，对其再进行滤波、重采样回原图大小；

步骤三，雾区校正；通过雾厚度图像HTM和所提取的fog_region二值图像，求取地物背景贡献值BACK，以此来修正HTM，达到保留正常区域的目的，最后基于辐射传输原理，在检测出的DCM雾区内，通过HTM去雾模型对原始影像进行校正处理。

2.根据权利要求1所述的适用于光学遥感影像的DCM-HTM去雾霾方法，其特征在于，在所述步骤一中，暗通道图像J^dark(x)的获取方法如公式所示；其中，J^c表示原始图像的每个通道，c表示R、G、B三通道，Ω(x)表示以像素x为中心的一个窗口；

定义J_i表示第i波段原始影像，DCM_i(x)为其对应的暗通道图像，其获取方法如公式 ${\mathrm{DCM}}_i\left(x\right)=\underset{y\&Element;\mathrm{\&Omega;}\left(x\right)}{\min }\left(J_i\left(y\right)\right)$ 所示；

定义总的概率密度为两个规则概率密度和乘积的结果，其中，x³是坐标空间部分，x^r是对应的颜色空间部分，h_s和h_r分别代表坐标空间和颜色空间的核函数带宽，c为归一化的相关系数，用于滤波并分割时所需设定的变量用h＝h(h_s，h_r，M)表示，其结果将包含像素点少于数值M的区域与它最相似的区域合并。

3.根据权利要求1所述的适用于光学遥感影像的DCM-HTM去雾霾方法，其特征在于，在所述步骤一中，经过mean-shift平滑处理后，利用公式 ${{\mathrm{DN}}^{\&prime;}}_{\mathrm{mean}\_\mathrm{shift}}=\left\{\begin{array}{cc}0,& {\mathrm{DN}}_{\mathrm{mean}\_\mathrm{shift}}<T\\ 1,& {\mathrm{DN}}_{\mathrm{mean}\_\mathrm{shift}}>T\end{array}\right.,$ 根据阈值T对暗通道图像二值化，将图像分为背景区域(DN′_{mean_shift}＝0)和雾区域(DN′_{mean_shift}＝1)；在[0，100]区间内直方图最大值i_histMax和最小值i_histMin处的DN值DN[i_histMax]和DN[i_histMin]来确定最佳阈值T，具体求值方法如公式T＝(DN[i_histMax]+DN[i_histMin])/2 DN[i_histMax]，DN[i_histMin]≠0所示。

4.根据权利要求1所述的适用于光学遥感影像的DCM-HTM去雾霾方法，其特征在于，在所述步骤二中，首先假定暗通道图像DCM宽高分别为W和H，根据所选取的不重叠窗口HTM_W(ω′，ω′)计算初始输出的HTM大小，即为DCM宽高的1/ω′；则从输出的HTM像元坐标点(x′，y′)出发，读取DCM起始点(ω′*x′，ω′*y′)，宽高均为ω′大小的窗口像元值，记为Ω(x′，y′)，最后取该窗口内的最小值作为当前点的HTM值，计算方法如公式HTM_i(x′，y′)＝min_y∈Ω(x′，y′){DCM_i(y)}所示；其中，DCM_i(y)为第i波段暗通道图像在窗口Ω(x′，y′)范围内的像元值，HTM_i(x′，y′)则表示对应的宽高为(W/ω′，H/ω′)的雾厚度图像。

5.根据权利要求4所述的适用于光学遥感影像的DCM-HTM去雾霾方法，其特征在于，对HTM进行滤波，选择中值滤波，模板大小为m×m，m∈[3，9]，来平滑矩阵窗口所产生的噪声并重采样成原图大小。

6.根据权利要求5所述的适用于光学遥感影像的DCM-HTM去雾霾方法，其特征在于，所述HTM包括三部分，分别是大气中气溶胶等的散射部分、一些高亮度地表以及雾的影响部分；这里将前两者归为大气的影响，以背景BACK表示，则修正后的雾厚度图像HIM′用 ${\mathrm{HTM}}^{\&prime;}=\left\{\begin{array}{cc}0,& \mathrm{HTM}\&le;\mathrm{BACK}\\ \mathrm{HTM}-\mathrm{BACK},& \mathrm{HTM}>\mathrm{BACK}\end{array}\right.$ 表示；则其中，为原始影像第i波段点(x，y)处的灰度值，为其校正的结果图像，HTM_i(x，y)为雾厚度图像。