CN103824260A - 一种高分辨率遥感图像薄雾快速去除技术 - Google Patents

Abstract

本发明针对波长在可见光与近红外范围的高分辨率遥感图像提供一种快速、自动的薄雾去除技术，该技术主要包括以下步骤：从目标图像本身估计薄雾的分布与厚度；在有与目标图像同一区域、获取季节相近且清晰的参考图像的情况下，利用参考图像更加精细地排除地表特性的影响，获取精细的薄雾数据；将估计出的薄雾数据代入薄雾成像模型，得到去雾后的图像。该算法利用二次大尺度中值滤波技术估算薄雾，处理流程简单、无需人工交互、处理速度快的特点，算法稳定，适用于薄雾厚度大、分布不均匀、不适合绝对辐射校正的高分辨率遥感数据的辐射处理。

Description

一种高分辨率遥感图像薄雾快速去除技术

技术领域

本发明涉及遥感图像薄雾去除技术，具体的说，涉及一种针对可见光到短波近红外波段(0.4～1μm)高分辨率遥感图像的薄雾快速去除技术。 

背景技术

遥感图像在成像过程中会受到大气分子、气溶胶和云粒子等大气成分的吸收与散射影响，为了获取地表真实地物的信息就要去除这些影响。目前遥感技术已经逐渐从定性发展到定量阶段，对于国外多数传感器数据，去除大气影响的首选方法是绝对大气校正(Absolute Atmospheric Correction)，绝对大气校正是基于大气辐射传输物理过程的，将遥感图像的DN(digital number)值转换为地表反射率或地表反射辐亮度的方法。研究者们发展了很多算法来消除太阳辐照度、大气和传感器差异带来的影响，例如：6S、MODTRAN等辐射传输模型。该类方法较复杂，需要利用一系列参数(例如，卫星过境时的地物反射率、大气能见度、太阳天顶角和卫星传感器标定参数等)。 

薄雾是大气影响严重到一定程度的表现，分布与厚度都不均匀的薄雾严重影响遥感图像的质量，对于目前很多高分辨率遥感图像，很难使用绝对大气校正来达到去除薄雾的目的，主要原因有2个，首先是绝对大气校正方法的限制，当薄雾厚度增加到一定程度(比如表现为气溶胶厚度超过1)，绝对大气校正依赖的6S、MODTRAN等辐射传输模型结果将产生很大误差。其次是高分辨率遥感图像的限制，这些数据波段设置往往只有可见光到短波近红外(0.4～1μm)，缺少绝对大气校正方法常用的长波近红外波段(2.1μm)，并且绝对大气校正所需要的参数获取困难。对于无法采用绝对大气校正的高分辨率遥感图像，目前多依赖基于图像图像本身的薄雾去除技术。商业软件ERDAS采用基于直方图匹配的方法，人工将图像划分为清晰和模糊区域，假设两类区域真实地表的直方图是相同的，匹配模糊区域与清晰区域的直方图达到薄雾去除的目的，方法虽然简单，但效果有限，并且需要人工交互。还有一些方法基于同态滤波、小波变换等，这类方法普遍存在处理效果有限、参数设置困难、运算速度较慢、很难适应大尺寸、大数据量遥感图像的去雾处理。 

目前中国大力发展卫星遥感技术，近年来发射了一系列的高分辨率遥感卫星，例如高分系列卫星、资源系列卫星与环境减灾系列卫星等，获取的高分辨率遥感图像数据量急速增加，随着近年来中国环境污染导致的雾霾现象急剧恶化，薄雾越来越多的出现在这些高分辨率遥感图像中，在无法采用绝对大气校正去除薄雾影响的情况下，研究一个有效的、自动的、快速的薄雾去除技术，满足海量高分辨率遥感图像去除薄雾的工程化要求成为目前国家的迫切需求。 

发明内容

本发明的目的是针对遥感应用，提供一种快速的高分辨率遥感图像薄雾去除技术，特别是对于那些可见光到短波近红外波段的高分辨率遥感图像的、厚度与分布都不均匀的薄雾的去除。本技术利用二次大尺 度中值滤波技术从图像本身估计薄雾的分布与厚度，并利用半物理的薄雾成像模型去除薄雾的影响。 

本发明的基本思路为：针对波长在可见光与近红外范围内的高分辨率遥感图像，先提取每个波段的暗通道数据，然后在暗通道数据上利用大尺度的中值滤波进行两次过滤，提取图像中具有大范围、高亮度、低频的信息来近似薄雾的分布与厚度，最后将估算出的薄雾数据代入半物理的薄雾成像模型，获得去除薄雾后的图像。本技术可以基于单张图像去除薄雾，估算出的薄雾有可能包含一些地物的信息，也可以利用相同区域的清晰图像作为参考数据，区分出地物的贡献，从而取得更好的去雾效果。 

所述的高分辨率遥感图像即包括卫星传感器数据，比如国产高分系列卫星、资源系列卫星与环境减灾系列卫星中的高分辨率多光谱传感器数据，也可以来自近地的航空传感器数据，比如飞机、无人机等飞行器平台，这些数据一般最多只有可见光到短波近红外的4个波段，定量化程度不高，不适宜进行绝对大气校正来去除大气影响。 

本发明的技术方案提供的高分辨率遥感图像薄雾去除方法，其特征在于包括以下实施步骤： 

A从目标图像本身估计薄雾的分布与厚度； 

B在有与目标图像同一区域、获取季节相近且清晰的参考图像的情况下，利用参考图像更加精细地排除地表特性的影响，获取更加精细的薄雾数据； 

C将估计出的薄雾数据代入薄雾成像模型，得到去雾后的图像； 

上述实施步骤的特征在于： 

步骤A中所述从目标图像本身估计薄雾的分布与厚度，主要依据的假设是薄雾在高分辨遥感图像中具有大面积、高亮度、低频的特性，利用该假设，通过图像处理的方法估算薄雾，具体是先计算图像每个波段的暗通道数据，然后基于暗通道数据，利用二次大尺度中值滤波技术获得薄雾的估计。步骤A中采用的二次大尺度中值滤波技术可以很好地适应高分辨率遥感图像的分布与厚度特征，且有计算量小、处理速度快的特点。大尺度中值滤波算法采用与滤波半径大小无关的快速算法。 

步骤B中所述的参考图像是与目标图像来自相同传感器，具有相同的地理区域，且获取季节相近，受大气影响小的、相对清晰的图像。由于步骤A中依据薄雾假设估算出的薄雾数据可能含有一定的地物信息，利用对参考图进行步骤A中的计算，获取地表估计，然后利用平滑窗口技术从步骤A中获得的薄雾估计数据中去除其中的地表信息，得到更加精细的薄雾估计。在没有参考图像的情况下，步骤B可以省略。 

步骤C中所述薄雾成像模型，来自定量遥感领域的大气辐射传输理论，描述了太阳光辐射与地表、大气作用后到达传感器的过程，该模型在步骤C中的计算不带物理单位，是一种半物理的模型。 

本发明与现有技术相比有如下特点：通过图像本身估算薄雾的分布与厚度，相比传统的基于同态滤波、小波变化等提取低频信息来去除薄雾的方法，该技术利用二次大尺度中值滤波技术获取图像中具有大面积、高亮度、低频信息作为薄雾的估计，具有处理流程简单、无需人工交互、运算速度快的特点，而且在算法的稳定性以及适用性上都具有优势。为中国目前那些定量化程度低、空间分辨率高的可见光到短波近 红外多光谱遥感图像薄雾去除提供了一种简单的解决方案。该技术结合图像分块处理技术及OpenMP加速技术后可完成大尺寸、海量遥感图像的工程化处理需求。 

附图说明：

图1是高分辨率遥感图像薄雾去除流程图 

图2是薄雾成像模型示意图 

图3是薄雾去除结果示意图，前三个波段(RGB)采用ENVI软件默认2％linear增强显示，左为某高分数据，右为泰国地球观测系统卫星数据 

图4是薄雾去除结果示意图，前三个波段(RGB)采用ENVI软件默认2％linear增强显示，左为资源3号卫星数据，右为环境减灾卫星A星CCD1数据 

具体实施方式：

本技术的思想是利用二次大尺度中值滤波技术估算薄雾的分布与厚度，然后代入薄雾成像模型获得去雾后的图像，其核心的薄雾假设主要有以下3点： 

1)图像中薄雾区域的亮度值大于清晰区域； 

2)图像中薄雾区域主要集中在低频区域，清晰地物主要集中在高频区域。 

3)图像中薄雾具有一定大小的面积分布，且分布不均匀，厚度从最厚处到清晰区域是逐渐减小的过程，没有突变边缘。 

采用本发明实现高分辨率遥感图像薄雾去除的实施例如图1所示，现结合附图对其进行描述。 

处理单元111计算暗通道。利用暗通道原理dark channel[He，CVPR2009]提取图像每个波段的低频信息。根据dark channel原理，在图像任一象素位置上，如果该位置的雾的浓度低，则对于不同波段之间在该位置的周围一定区域内(比如5×5象素区域内)往往可以找到一个位置具有较低的象素值，即所谓的该位置的暗通道值，由暗通道值组成的新的波段可以很好地描述包括薄雾在内的低频数据的分布及厚度。记W(x，y)为遥感图像某个波段在像素位置(x，y)处一定局部范围内的最小值组成的新的波段数据，即该波段的暗通道。 

处理单元112二次大尺度中值滤波。对W(x，y)进行进一步处理，因为W(x，y)的低频信息还包含大量具有低频特性的地物信息，我们需要分离出更大尺度的低频信息做为薄雾的估计。我们采用基于大滤波半径的中值滤波方法，公式如下： 

A(x，y)＝median_r(W)(x，y) 

B(x，y)＝p(A(x，y)-median_r(|W-A|)(x，y)) 

这里median表示中值滤波，r为中值滤波窗口的半径，根据薄雾在高分辨率遥感图像中的尺度，r值往往取值较大，一般从20象素到100象素都有可能。上述第一个公式对W(x，y)进行大尺度核的中值滤波 后，得到的是得到A(x，y)为局部平滑的图像，但又保持了部分的跃迁边缘，而减去二次平滑项median_r(|W-A|)(x，y)后。而上述第二个公式便主要是具有厚度渐变特性的薄雾了，公式中的p为取值0～1的乘子，取值越大，估算的V(x，y)值越大，表示雾霾越厚，通过实验发现p取值一般为0.8～0.95，我们算法中默认取值0.95。 

最后，考虑到大气光的取值范围不会超过暗通道数据，于是有0＜V(x，y)＜W(x，y)，最终： 

V(x，y)＝max(min(B(x，y)，W(x，y))，0) 

估计出的V(x，y)便可以作为具有大面积、高亮度、低频特性的薄雾。 

一般图像处理中的中值滤波运算速度是非常快的，这主要是因为通常采用的滤波窗口都很小，比如3×3象素、6×6象素，但是随着滤波窗口半径r的增大，对窗口内所有象素的排序运算量也增加，中值滤波的计算量快速增加。采用大尺度的中值滤波提取图像中的低频信息，就需要一种运算时间不随滤波窗口变大而增加的、尺度无关的中值滤波算法。同时随着滤波窗口的增大，正方形窗口已经不适用，需要采用尽量接近圆形的滤波窗口。采用Perreault的中值滤波算法，该算法不论滤波窗口半径r如何变化，算法复杂度始终为O(1)，并且采用圆形滤波窗口。Perreault提供了算法的标准C开源代码，网址：http://nomis80.org/ctmf.html。但是该开源算法只支持8位Byte型数据处理，Perreault利用8位数据生成的存储直方图分级为256。我们需要中值滤波算法支持16位unsigned short型数据，因为目前很多高分辨率遥感图像的存储数据类型是Uint16。根据Perreault的方法，16位数据的直方图需要65536分级，这会导致算法速度降低256倍，我们采用2ⁿ倍数减小16位直方图分级，提高中值滤波计算速度。在Perreault源码的基础上，我们添加了对16位Uint16数据的支持，并且可以设置直方图分级数据。 

处理单元113地表信息分离。在有参考图像的情况下，对参考图像进行和处理单元111、处理单元112相同的计算，得到地表估计数据，记V_L(x，y)，记通过目标图像得到的薄雾估计数据为V_H1(x，y)，为了保证V_l(x，y)与V_H1(x，y)在相同的数量级上，统计V_H1(x，y)的直方图，然后以V_H1(x，y)数据前10％的像素所在位置，主要是图像中清晰像素的位置，找到对应V_L(x，y)位置上的像素点，然后利用这两组像素值利用最小二乘原理求出V_L(x，y)到V_H1(x，y)的线性变换关系，利用该线性变化修正V_L(x，y)的值。 

地表信息分离后得到的薄雾估计V_H2(x，y)按照以下公式计算： 

$V_{H2}(x,y)=\left\{\begin{array}{cc}{V}_{H1}(x,y)-V_L(x,y),& \mathrm{if}{V}_{H1}(x,y)\≥V_L(x,y)\\ 0,& \mathrm{if}{V}_{H1}(x,y)\≤V_L(x,y)\end{array}\right.$

按照上述公式求出V_H2(x，y)，为了保证数据的平滑，采用5×5的平滑滤波窗口对V_H2(x，y)进行一 次平滑处理，得到最终的薄雾估计。 

处理单元114基于薄雾成像模型的薄雾去除。薄雾成像模型可以认为是辐射传输方程的一种简化情况，该模型将传感器接收到的辐射分成两部分，一部分为太阳光经过大气衰减后到达地面，经地面反射后又经过一次大气的衰减后到达传感器的辐射亮度。另外一部分为太阳光照射到大气分子上反射后到达传感器的辐射亮度。对于第一部分，为与地表作用的那部分太阳辐射。假定太阳辐射强度为I(x，y)，则经过指数衰减到达地面的辐射强度为I(x，y)e^-kλU，假定在地表(x，y)位置的二向反射分布函数为b(x，y)，则经过地表反射后辐射亮度变为I(x，y)e^-kλUb(x，y)，再经过大气衰减后最终到达传感器的辐射亮度为I(x，y)e^-kλub(x，y)e^-kλU。对于第二部分，为与大气作用的那部分太阳辐射。太阳光I(x，y)经过大气作用后衰减为I(x，y)e^-kλU，最终到达传感器的辐射强度为I(x，y)(1-e^-kλU)。将两部分加起来便是传感器最终接收到的辐射L(x，y)。 

L(x，y)＝I(x，y)e^-kλub(x，y)e^-kλU+I(x，y)(1-e^-kλU) 

上述公式是辐射传输方程的一种简化形式，但对于定量程度不高的高分遥感数据，求解还是太困难了，需要简化，于是记J(x，y)＝I(x，y)e^-kλUb(x，y)，则J(x，y)为地表的真实辐射亮度，也就是我们希望提取的去雾后图像。于是： 

L(x，y)＝J(x，y)e^-kλu+I(x，y)(1-e^-kλU) 

该公式示例见图2，公式中描述大气作用的指数项e^-kλU求解比较困难，因为对于定量程度不高的高分遥感数据，参数k、λ、U很难给出相对精确的值，只能通过图像的特征去估计，而通过图像特征去估计e^-kλU本身是困难的，这里不考虑直接去估算e^-kλU的值，而是估算去太阳辐射大气分子上反射后到达传感器的辐射亮度，于是记V(x，y)＝I(x，y)(1-e^-kλU)，此时公式变为： 

$L(x,y)=J(x,y)(1-\frac{V(x,y)}{I(x,y)})+V(x,y)$

故此时地表真实的辐射亮度J(x，y)表示为： 

$J(x,y)=\frac{L(x,y)-V(x,y)}{1-\frac{V(x,y)}{I(x,y)}}$

这样J(x，y)的求解就只含有对V(x，y)与I(x，y)的估计。V(x，y)便是估算出的薄雾数据，I(x，y)对于高分遥感数据的每个波段近似固定值，一般由传感器生产方提供，即使没有提供，也可以用每个波段亮 度的最大值近似代替。 

薄雾成像模型虽然来源于大气辐射传输理论，但计算过程很难量化到绝对物理量的计算，是在图像数据的象素DN值上的计算，没有物理单位，相对于定量遥感中严格的辐射传输物理模型，将这里的薄雾成像模型称为半物理模型。 

将估算出的薄雾数据V(x，y)代入薄雾成像模型，得到去雾后的图像。 

本发明实现的高分辨率遥感图像薄雾去除算法计算简单，结合尺度无关中值滤波技术后，处理速度远优于基于同态滤波、小波变换等薄雾去除算法，且无需复杂参数设置，针对具体卫星传感器的图像处理可以实现完全自动化的批处理。与基于绝对大气校正的大气作用消除做法不同，本方法定量化程度低，计算数据不带物理单位，数据在去除薄雾的同时，地物的光谱特性会有一定程度的改变，因此本方法处理后的数据不适宜遥感定量应用，适宜遥感制图、分类等应用。 

率的处理流程比传统的绝对辐射校正流程要简化。传统的绝对辐射校正流程，需要辐射定标与绝对大气校正两步，并且这两步的处理需要大量的参数，包括定标系数、传感器参数与大气参数等，并且基于辐射传输理论的绝对大气校正运算复杂，处理速度慢。本发明的处理流程，处理流程只有相对辐射校正一步，输入的参数也少，无需定标参数及精确的大气参数。并且本发明的计算过程简单，处理速度快。 

本发明对四种高分辨率遥感图像的适用性分析如下： 

为了分析算法的适用性，列举四种不同空间分辨率的遥感图像数据的薄雾去除结果，分别为：空间分辨率8米的中国某高分系列卫星(GF)、空间分辨率2.5米的泰国地球观测系统卫星(THEOS)融合数据、空间分辨率6米的资源三号卫星(ZY-3)，与空间分辨率30的环境减灾卫星(HJ-1A CCD)。具体效果见图3与图4。图3左为GF数据，图3右为THEOS数据，图4左为ZY-3数据，图4右为HJ-1A CCD数据。图中数据波段按照RGB彩色显示，显示方式本身经过2％线性拉伸增强。几张遥感图像中薄雾的类型、分布与厚度存在差异，算法都可以在很大程度上去除薄雾的影响。 

本发明的一个实例在PC平台上实现，支持Windows与Linux操作系统，目前支持大多数高分辨率遥感图像薄雾去除处理，支持Byte型与Uint16型数据处理，唯一的参数是中值滤波半径，该参数在确定了具体传感器后，根据图像的空间分辨率便可以基本确定数值，用户可以在开发的软件中的配置文件中添加不同传感器的配置信息，也可根据具体数据手动更改。薄雾去除效果受数据本身影响，在通过大量的数据处理与验证的基础上，确定了常用高分辨遥感图像对应的中值滤波半径，对这些数据的处理做到了自动化。算法计算原理简单，计算量小，处理速度快，结合简单的OpenMP加速技术后，在普通单机32位、Intel core i7-3770，CPU3.4GHz配置下，对于10000×10000像素尺寸4波段图像，数据类型为Byte型时处理速度接近数据IO数据，对于Uint16数据类型，平均处理时间小于5分钟，并且算法采用了图像分块处理技术，支持超大尺寸遥感图像处理，考虑到大尺寸中值滤波对边缘象素处理的差异，图像分块时，相邻块保 持一个中值滤波窗口长度的重叠区域，便保证分块处理的结果同整体处理的完全一致性。软件鲁棒性较强，满足高分辨率遥感图像去雾自动、快速处理的工程化需求。 

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。 

Claims (5)

1.一种高分辨率遥感图像薄雾快速去除技术，该方法针对波长在可见光到短波近红外范围内的高分辨率遥感图像中存在的厚度与分布都不均匀的薄雾的去除，其特征在于：

首先从目标图像本身估计薄雾的分布与厚度，然后在有与目标图像同一区域、获取季节相近且清晰的参考图像的情况下，利用参考图像更加精细地排除地表特性的影响，获取更加精细的薄雾数据，最后将估计出的薄雾数据代入薄雾成像模型，得到去雾后的图像。

2.根据权利要求1中所述的技术方案，其特征在于：

从目标图像本身估计薄雾的分布与厚度，具体流程是先计算图像每个波段的暗通道数据，然后基于暗通道数据，利用二次大尺度中值滤波技术获得薄雾的估计，对暗通道数据计算方法不做权利要求。

3.根据权利要求1中所述的技术方案，其特征在于：

利用二次大尺度中值滤波技术获得薄雾的估计，其核心公式为：

A(x，y)＝median_r(W)(x，y)

B(x，y)＝p(A(x，y)-median_r(|W-A|)(x，y))

V(x，y)＝max(min(B(x，y)，W(x，y))，0)

其中median表示中值滤波，r为中值滤波窗口的半径，取值较大，一般20～100象素，W(x，y)为每个波段的暗通道数据，即遥感图像某个波段在像素位置(x，y)处一定局部范围内的最小值组成的新的波段数据，p为取值0～1的乘子，用于描述薄雾的厚度，V(x，y)为估算出的薄雾数据。

4.根据权利要求1中所述的技术方案，其特征在于：

利用二次大尺度中值滤波技术获得薄雾的估计，采用支持16位unsigned short型数据的中值滤波算法，算法复杂度与滤波半径无关，终为O(1)，并且采用圆形滤波窗口，对支持8位unsigned char型数据的中值滤波算法不做权利要求。

5.根据权利要求1中所述的技术方案，其特征在于：

将估计出的薄雾数据代入薄雾成像模型，得到去雾后的图像，薄雾成像模型公式为：

$J(x,y)=\frac{L(x,y)-V(x,y)}{1-\frac{V(x,y)}{I(x,y)}}$

其中J(x，y)为去雾后图像，V(x，y)为估算出的薄雾数据，I(x，y)对于高分遥感数据的每个波段近似固定值，一般由传感器生产方提供，没有则以用每个波段亮度的最大值近似代替。

Images