1 引言

雾霾由空气中的水汽、灰尘、烟、雾等微小颗粒组成,它会影响大气透过率,导致大气的能见度下降^[1]。遥感卫星成像时,雾霾的影响近似于在遥感卫星成像通道中加入了一层毛玻璃,造成雾状模糊,导致图像对比度下降、清晰度降低,从而影响图像解译和进一步应用^[2-3]。研究雾霾(包括均匀薄云)对遥感成像的影响机制,可减小或消除这类大气成分对遥感图像质量的影响,是提高遥感图像质量和应用价值的重要途径。

遥感卫星成像通道中雾霾的影响机制不同于近景摄影测量中地面拍摄时的雾霾,遥感中雾霾垂直于大气方向,而地面拍摄时雾霾平行于大气方向^[4],因此处理方法也不尽相同。近景摄影测量中常用的图像去雾霾方法可分为两大类:1)采用图像增强的方法,主要有暗通道先验理论、同态滤波及增强同态滤波、Retinex等方法,也有学者将上述方法用于遥感图像的去雾霾处理,但这些方法存在图像信息损失或过度增强等现象,影响遥感图像的解析^[5-6];2)采用图像复原的方法,以雾霾天降质图像的物理模型为基础,根据一定的假设或先验知识对雾霾图像进行复原^[7-8],这种方法的局限性在于雾霾降质图像的物理模型精度难以适应复杂遥感成像的需求。

目前遥感图像的去雾霾处理方法多数基于暗像元或暗通道相关方法,如:Makarau等^[9-10]提出利用暗像元来去除雾霾对图像DN(Digital Number)值影响的方法,并将该方法成功应用于Worldview-2和Landsat-8 OLI等卫星图像;Scaramuzza等^[11]基于Landsat-8 OLI图像的研究提出了AT-ACCA(Artificial Thermal-Automated Cloud-Cover Assessment)方法;许继伟等^[12]提出在暗通道的基础上构建雾霾图谱的方法,并利用Landsat-8 OLI图像进行验证。国产高分辨率卫星(如GF-2、Superview-1等)的谱段设置一般为红、绿、蓝和近红外4个多光谱谱段和1个全色谱段^[13],难以比拟Worldview-2和Landsat-8 OLI的8个多光谱谱段和1个全色谱段,上述方法直接移植在国产高分辨率遥感图像去雾霾处理时难免存在局限性。

本文基于常用的暗通道处理原理,在Makarau方法的基础上,提出了一种充分利用国产高分辨率遥感卫星谱段特性的去雾霾方法,并通过GF-2和Superview-1图像进行验证。

2 算法原理

基于图像暗通道原理,将Makarau方法扩展到遥感图像去雾霾的处理。遥感图像受大气的影响可分为两部分:第一部分为大气成分中气溶胶、水汽等对辐射的散射^[14]、吸收以及程辐射的影响,其处理方法为基于辐射传输模型的大气校正;第二部分为雾霾(包括薄云)的影响,其处理方法是基于暗通道的去雾霾处理^[15]。根据上述介绍,遥感图像DN值组成可表示为

Dtotal=Dtarget+Dareo+wat+Dhaze,(1)

式中:D_total为图像总DN值;D_target为目标真实的DN值;D_areo+wat为气溶胶和水汽等影响的DN值,通过大气校正去除;D_haze为雾霾影响的DN值,通过所提方法去除。

所提去雾霾方法的基本思想可描述为:通过暗通道分析雾霾对图像DN值的影响,从原始图像DN值中去除这一部分雾霾DN值的影响,恢复被摄目标景物的真实DN值信息。这种直接从原始图像DN值中扣减雾霾DN值的方法类似减去背景DN值处理,由于大气成分之间多为耦合作用,D_haze和D_areo+wat之间并非完全独立,扣减的这部分D_haze可能包含部分气溶胶和水汽的影响,尤其在可见近红外谱段主要包含气溶胶光学厚度的影响。因此,需要进行气溶胶补偿(下文将对大气成分的补偿称为气溶胶补偿),避免过度扣减背景而造成目标信息丢失。

3 处理流程

3.1 图像分割

遥感图像中景物种类繁多且特征不一,灰度图像中景物亮度特征都表现为DN值,很难直接区分背景地物、自身亮度较高的地物(此类地物一般地表反射率较大,图像DN值较高,如屋顶、挡风玻璃、沙地等,下文中称之为亮目标)和亮度偏高的雾霾区域。为了精确定位雾霾区域,对图像进行分类。分类的目的主要是将容易与雾霾相混淆的亮目标进行分类处理,避免由误判导致图像失真。

图像分类的方法有多种,所提方法采用最常用的mean-shift算法^[16],定义其算核M(x)为

M(x)=∑k=1nxkgk(x)∑k=1ngk(x)-x,(2)

式中:n为采样点数;k为采样点序号;x_k为一个L维空间中n个采样点中的第k点;g_k(x)为核函数(通常为高斯核函数);x为某个样本点。其中,高斯核函数g_k(x)可表示为

gk(x)=exp-‖xsk‖2hsexp-‖xrk‖2hr,(3)

式中:上标s和r分别表示采样区间和灰度; xsk为坐标; xrk为灰度信息;h_s和h_r分别为采样区间大小和灰度信息偏差,称为mean-shift分割算法的带宽,带宽过小易造成欠分割,带宽过大易造成过分割^[17]。结合遥感图像景物种类多且变化复杂的特点,对mean-shift算法进行改进,将其核函数中通常固定的两个带宽量改进为自适应的带宽量,即带宽改进为分类区域大小及灰度的函数,其表达式为

hs=1MN∑a=1MNak(a),(4)hr=1MN∑i=1M∑j=1Ngra[f(i,j)],(5)

式中:k(a)为分类区域某一个样本点的归一化距离;a为对应采样点的位置序号;M、N分别为图像长和宽的像素数;f(i,j)为图像(i, j)点对应的灰度值;gra[·]表示求灰度梯度运算。

图像分割的主要目的是将亮目标分割出来,对mean-shift分割后的图像用阈值进行分类,将分割图像的均值作为阈值、大于阈值的像元部分作为亮目标,即

Dbr_obj(i,j)=0,f'(i,j)≤T1,f'(i,j)>T,(6)

式中:D_{br_obj}(i,j)表示亮目标图像的DN值;f'(i,j)为分割后图像中(i, j)像素的灰度值;T为亮目标判断阈值,其大小等于分割后图像DN值的均值。

3.2 暗通道图像处理

暗通道图像处理包括两个步骤:1)构建暗通道检测掩模,将暗通道与非暗通道进行区分;2)采用暗通道窗口检测出雾霾图像。

步骤2)主要包含两个部分。

(1) 检测暗通道图像。根据Chavez的研究^[15],遥感图像中反射率低的裸露土壤、浓密植被暗通道像元等更易受到雾霾的影响,因此暗通道检测掩模可采用不重叠的求最小值掩模,对图像进行逐窗口处理,每一个窗口均输出一个最小值。Makarau等^[9-10]的研究表明窗口大小为对应遥感图像空间分辨率的15~20倍,通过3次插值将图像恢复至原始图像大小,采用中值滤波的方法抑制插值后可能出现的振铃现象。通过阈值计算,将小于经过掩模检测的图像均值的像素部分作为暗通道像素:

f″(i,j)=Wmin*f(i,j),(7)Cdark=0,f″(i,j)>Tdark1,f″(i,j)≤Tdark,(8)Cno_dark=1,f″(i,j)>Tdark0,f″(i,j)≤Tdark,(9)

式中:*为卷积运算;f″(i,j)为掩模卷积后的图像灰度值;W_min为计算最小值的掩模窗口;C_dark为二值化的暗通道图像灰度值;C_{no_dark} 为二值化的非暗通道图像灰度值;T_dark为暗通道阈值,其大小为f″(i,j)的均值。暗通道图像主要用于雾霾校正,而非雾霾图像则主要用于气溶胶等大气成分贡献值的补偿。

(2) 构建雾霾图像。与暗通道检测方法类似,利用更小的求最小值掩模对原始图像进行逐窗口处理,每个窗口均输出对应最小值,并将此最小值作为该区域的雾霾对DN值的贡献。Makarau等^[9]依据对30 m分辨率的Landsat-8 OLI卫星图像的研究提出窗口大小为遥感图像空间分辨率的3~5倍时可达到最佳检测效果。考虑到实际中一幅遥感图像景物可能出现较大的高亮物体或区域,尤其对亚米级空间分辨率的遥感图像,采用Makarau等^[9]提出的最大5倍窗口的掩模,也会出现窗口区域内景物平均亮度很高的情况,使得掩模输出的最小值无法作为该区域内的雾霾DN值,这种固定掩模窗口大小的方法应用于高分辨率遥感图像时存在很大的局限性。因此,提出了一种自适应雾霾区域检测掩模的方法,以3 pixel×3 pixel作为初值掩模窗口大小,每一次进行掩模卷积时均对输出的最小值进行判断后再输出。当该值小于原始图像灰度均值时,将其作为有效雾霾DN值输出;当该值大于原始图像灰度均值时,此时检测窗口自动增大1倍并继续检测,直到满足条件为止。雾霾DN值可表示为

Dhaze=Wminnew*f(i,j),(10)

式中: Wminnew为初始大小为3 pixel×3 pixel的求最小值掩模窗口,上标new表示新掩模窗口。

3.3 雾霾图像补偿与校正

第2节中已经介绍,如果直接从原始图像DN值中减去由暗通道处理得到的雾霾DN值,将会造成过校正,从而导致图像失真。出现该现象的主要原因有两个:1)通过暗通道处理得到的雾霾图像可能包含部分亮目标的影响;2)由于大气成分之间的耦合作用,雾霾DN值中包含气溶胶等贡献。因此还需对雾霾图像进行亮目标校正并对气溶胶等大气影响进行补偿。

亮目标补偿的原理是通过比对第3.1节中得到的亮目标图像与第3.2节中的雾霾图像,对二者重叠区域进行亮目标补偿。重叠区域包含了亮目标和雾霾的共同贡献,对遥感图像进行处理时一般采用插值的方法对异常像素进行处理(此处亮目标即可理解为异常像素),常用的插值方法包括邻近像元插值法、双线性内插法以及邻域均值插值法等,一般可用作一个或几个像素的插值。本研究中高分辨率遥感图像分割得到的亮目标可能占数十甚至上百个像素,采用常用的插值方法可能导致无法对亮目标进行有效补偿,尤其是亮目标内部的像素,其DN值依旧偏高,常用的插值方法易导致校正后的高分辨率图像出现偏色。因此提出一种三角形插值的方法,利用重叠区外3个临近非亮目标的像素DN值,对亮目标区域进行逐像素插值,其原理示意图如图1所示。

图 1. 三角形插值示意图

Fig. 1. Schematic of triangulation interpolation

下载图片 查看所有图片

图1中黑色像素为雾霾图像中的亮目标,假定亮目标像素1的坐标为(i, j),亮目标像素2的坐标为(i+1, j+3),则像素1、2的DN值计算公式分别为

D1=D(i,j-1)+D(i-1,j)+D(i,j+1)3,(11)D2=D(i+1,j+2)+D(i+2,j+3)+D(i-1,j+4)3。(12)

对气溶胶补偿的原理是基于非雾霾区域主要受气溶胶影响而不受雾霾影响,利用第3.2节中得到的C_{no_dark}进行补偿,计算公式为

D'hz_img=Dhz_img-|(Dno_hz_img-Dhz_imgno_hz_img)|,(13)

式中:D'_{hz_img}为气溶胶补偿后的雾霾图像DN值;D^no_hz_img 为原始图像中非雾霾区域DN值(由第3.2节中C_{no_dark}得到); Dhz_imgno_hz_img为雾霾图像中非雾霾区域的DN值。

3.4 结合国产高分辨率卫星谱段特点的处理方法

目前国产高分辨率卫星的谱段设置多为4+1,即4个可见近红外多光谱谱段(蓝、绿、红、近红外)和1个全色谱段。研究表明^[18],在常用的可见近红外谱段(450~900 nm)中,较短波长的探测谱段更易受到大气吸收和散射的影响,大气透过率较低,随波长的增加,大气透过率逐渐增加。换言之,光学遥感中一般谱段波长越长,受气溶胶和雾霾等大气成分的影响越小,如近红外和红谱段,而波长较短的蓝谱段则更容易受到这些大气成分的影响。

据此,针对国产高分辨率遥感卫星的谱段特点,利用受大气影响更大的蓝谱段进行全局雾霾图像的检测,采用同样的方法得到其他谱段的雾霾图像。由于敏感性不同,采用(13)式对蓝谱段外其他谱段进行去雾霾处理时容易造成过校正,导致图像信息丢失和色彩失真。因此,利用线性拟合的方法拟合出各谱段雾霾图像与蓝谱段雾霾图像之间的关系,从而推出其他多光谱谱段的雾霾图像,并推广到全色谱段。对应I谱段的拟合系数可表示为

k=linear_fitting(DIhaze,Dbhaze),(14)

式中:I为谱段号;k为对应I谱段的拟合系数; DIhaze为对应I谱段的雾霾图像; Dbhaze为蓝谱段雾霾图像;linear_fitting(·)表示线性拟合运算。通过迭代的方法,得到各谱段稳定的k值。

利用受大气影响更小的红谱段进行图像分割和亮目标提取(此处采用红谱段而不采用波长更长的近红外谱段的依据是近红外谱段除了反射辐射外还存在部分热辐射,可能会带来温度的影响,而导致对亮目标的误判),完成亮目标校正。

4 实验验证与分析

实验对象为两颗国产亚米级高分辨率遥感卫星图像,即GF-2和Superview-1图像。两颗卫星参数如表1所示。其中,GSD为地面采样距离,即空间分辨率;PAN表示全色谱段;B1~B4表示4个不同波段;MS表示多光谱谱段。从表1中可知其谱段设置为典型的大气窗口,可避免臭氧、水汽等大气成分的额外影响,对图像质量产生影响的大气成分主要是雾霾粒子和气溶胶分子。

表1中GF-2图像为2014-10-06对深圳成像的图像,Superview-1图像为2018-03-25对上海成像的图像,原始图像为经过相对辐射校正预处理后的一级图像产品。由于一幅遥感图像通常达到数亿像素,进行算法验证实验时在原图中选取部分细节图,如图2所示。

原图中,由于雾霾的影响,通过主观观测可发现毛玻璃背景的能见度不高。出现这种现象的原因是雾霾使得图像产生了一定程度的模糊。采用所提方法分别对GF-2和Superview-1图像进行处理,处理过程中关键图像包括的短波长谱段的雾霾图像(蓝谱段)、长波长谱段的亮目标图像(红谱段)、暗通道图像和非暗通道图像分别如图3、4所示。

图 2. 原始图像。(a) GF-2 全色图像;(b) GF-2 多光谱图像;(c) Superview-1 全色图像; (d) Superview-1多光谱图像

Fig. 2. Original images. (a) GF-2 panchromatic image; (b) GF-2 multispectral image; (c) Superview-1 panchromatic image; (d) Superview-1 multispectral image

下载图片 查看所有图片

图 3. 对GF-2进行去雾霾处理的图像。(a)蓝谱段雾霾图像;(b)红谱段亮目标图像;(c)暗通道图像;(d)非暗通道图像

Fig. 3. GF-2 images after dehazing process. (a) Haze image in blue spectral band; (b) image of bright object in red spectral band; (c) dark channel image; (d) non-dark channel image

下载图片 查看所有图片

图 4. 对Superview-1进行去雾霾处理的图像。(a)蓝谱段雾霾图像;(b)红谱段亮目标图像;(c)暗通道图像;(d)非暗通道图像

Fig. 4. Superview-1 images after dehazing process. (a) Haze image in blue spectral band; (b) image of bright object in red spectral band; (c) dark channel image; (d) non-dark channel image

下载图片 查看所有图片

采用所提方法和Makarau等^[9]提出的方法同时对图像进行处理,GF-2和Superview-1图像处理前后对比图如图5、6所示。

对所提方法和Makarau方法处理前后的图像质量进行评价。雾霾的作用机制相当于增大了大气点扩散函数及大气吸收,结果体现在图像中景物边缘宽度增加,边缘锐利程度降低,图像清晰度下降。因此,采用图像平均灰度梯度、图像边缘能量、图像对比度、图像锐度这几个基于景物边缘信息变化的典型高分辨率遥感图像质量参数进行评价。为方便对比,每组图像质量参数计算结果均以该组中最大值为分母进行归一化处理,结果如表2~5所示。

为了更形象地研究去雾霾处理校正前后的遥感图像质量变化,将GF-2和Superview-1卫星图像每一种质量参数各谱段的结果进行对比,如图7、8所示。

图 5. GF-2去雾霾处理前后的图像。(a)原始全色图像;(b) Makarau方法处理后的全色图像;(c)所提方法处理后的全色图像;(d)原始多光谱真彩色图像;(e) Makarau方法处理后的真彩色图像;(f)所提方法处理后的真彩色图像

Fig. 5. GF-2 images before and after dehazing processing. (a) Original panchromatic image; (b) panchromatic image after processing by Makarau's method; (c) panchromatic image after processing by proposed method; (d) original multispectral true color image; (e) true color image after processing by Makarau's method; (f) true color image after processing by proposed method

下载图片 查看所有图片

图 6. Superview-1去雾霾处理前后的图像。(a)原始全色图像;(b) Makarau方法处理后的全色图像;(c)所提方法处理后的全色图像;(d)原始多光谱真彩色图像;(e) Makarau方法处理后的真彩色图像;(f)所提方法处理后的真彩色图像

Fig. 6. Superview-1 images before and after dehazing processing. (a) Original panchromatic image; (b) panchromatic image after processing by Makarau's method; (c) panchromatic image after processing by proposed method; (d) original multispectral true color image; (e) true color image after processing by Makarau's method; (f) true color image after processing by proposed method

下载图片 查看所有图片

图 7. 去雾霾处理前后GF-2图像质量参数的对比。(a)平均灰度梯度;(b)边缘能量;(c)对比度;(d) 锐度

Fig. 7. Comparison histogram of GF-2 image quality parameters before and after dehazing processing. (a) Average gray-scale gradient; (b) edge energy; (c) contrast; (d) sharpness

下载图片 查看所有图片

图 8. 去雾霾处理前后Superview-1图像质量参数的对比。(a)平均灰度梯度;(b)边缘能量;(c)对比度;(d)锐度

Fig. 8. Comparison histogram of Superview-1 image quality parameters before and after dehazing processing. (a) Average gray-scale gradient; (b) edge energy; (c) contrast; (d) sharpness

下载图片 查看所有图片

由表2~5以及图7、8可知:利用所提方法对遥感图像进行去雾霾处理后,图像质量有了明显的提升;相比Makarau等^[9]提出的方法,所提方法可得到更好的校正效果。

一般地,去雾霾过程中亮目标易与雾霾区域发生误判而导致光谱色彩特征失真,下文对Makarau等^[9]提出的方法和改进分割与亮目标校正后的所提方法的处理效果进行对比。部分亮目标局部示意图如图9所示。

图 9. 去雾霾处理前后局部亮目标色彩示意图。(a)原始图像;(b) Makarau方法处理后的图像;(c)所提方法处理后的图像

Fig. 9. Schematics of color of local bright object before and after dehazing processing. (a) Original image; (b) image processed by Makarau's method; (c) image processed by proposed method

下载图片 查看所有图片

图9中红色虚线框中表示发生偏色现象的地物。通过对比图像中局部高亮目标可知,原始图像中高亮度像素区域在Makarau方法处理后发生偏色,而所提方法处理后的图像很好地保持了原图色彩光谱特征。所提对图像亮目标分割及补偿的方法对图像细节和色彩的保持起到了很好的作用。

上述实验对象均为不同卫星一景图像中局部区域的特征目标,下图为对全景图像进行处理的结果,以GF-2卫星2014-10-06得到的图像为例。

所提算法应用于全景图像处理时依旧显示出了很好的应用效果,很好地去除了图像中的雾霾。不足之处在于:1)对大面积水体、浓密植被等低反射率目标区域的处理效果欠佳,出现了偏色现象;2)基于Matlab平台计算时计算效率不理想,对上述7300 pixel×7000 pixel的多光谱图像的处理总耗时约为3.6 h(上文中开窗处理局部图像时,每个窗口区域图像平均耗时为十几秒到数十秒)。

所提方法主要应用于谱段设置为红(0.63~0.69 μm)、绿(0.52~0.59 μm)、蓝(0.45~0.52 μm)及近红外(0.77~0.89 μm)的遥感图像,适用于国内大部分高分辨率系列卫星和资源系列卫星,如典型GF-1及其后续星、ZY-3(01/02)星,以及中等分辨率的HJ-1卫星等。实验结果表明:1)当雾霾严重到遮挡地物时,雾霾的影响几乎和云一样,采用所提方法很难较好地恢复地物细节及光谱信息;2)在全景图像去雾霾实验中,在低反射率区域(如实验中提到的浓密植被、水体等)存在偏色现象,原因在于不同谱段对浓密植被、水体等地物的响应率差异较大,尤其是红、蓝谱段响应的差异,在去雾霾处理后合成真彩色图像易出现偏色现象,低反射率区域的偏色基本表现为偏红。全景图像的去雾霾处理效果如图10所示。

图 10. 全景图像的去雾霾处理。(a)原始图像及细节;(b)去雾霾处理后的图像及细节

Fig. 10. Dehazing processing of whole scene image. (a) Original image and some details; (b) dehazing image and some details

下载图片 查看所有图片

5 结论

针对国产高分辨率遥感图像受雾霾影响发生图像质量退化的问题,在暗通道图像去雾方法的基础上提出了一种结合国产遥感卫星谱段特性的遥感图像去雾霾方法。以GF-2和Superview-1图像为实验对象,对所提方法的有效性进行了验证。实验结果表明:所提方法很好地去除了由雾霾造成的图像模糊;相比Makarau提出的去雾霾方法,所提方法对图像质量的提高程度更大,并且校正后的图像很好地保持了原始目标光谱色彩特性。大气校正是实现定量遥感的重要途径,随着高分辨率遥感数据的定量化应用需求越来越大,去雾霾处理后的遥感图像可作为大气校正的输入,从而提高高分辨率遥感图像的应用能力。所提方法对不规则分布的薄云及较严重的雾霾校正效果不佳,且数据量较大时影响计算效率,后续将针对薄云去除以及运算效率的提高进行进一步研究,为后期工程化应用在国产高分辨率影像预处理环节中提供一定的技术基础。