汪月云,黄 微,王 睿

(上海大学 通信与信息工程学院，上海 200444)(*通信作者电子邮箱yueyunW@163.com)

0 引言

光学遥感影像往往受到薄云污染的影响，造成影像地物细节对比度下降，这给影像的应用[1-3]带来了很多困难。因此，如何有效去除降质影像的薄云影响从而提高影像的利用率是遥感影像应用的一项重要研究。

目前，光学遥感影像薄云去除主要有频域法和空域法两类。频域法去云主要根据薄云信息处于影像低频区域的原理，所以可以设计合适的高通滤波器以滤除低频成分从而有效去云。比如:阎庆等[4]结合小波变换和同态滤波、冯春等[5]利用同态滤波都达到了有效去云效果。

空域法去云将薄云污染分为乘性污染和加性污染两种模型。乘性污染经典的去云方法有暗通道先验法[6-7]和薄云优化变换(Haze Optimized Transformation, HOT)法[8-9]。暗通道先验法通过对清晰影像库统计得到暗通道分布规律，即清晰影像至少在某一个颜色通道中有灰度值接近于0的像素，利用该先验规律得出云模型的大气光和透射率从而进行薄云去除。Pan等[10]研究发现无云遥感影像的暗通道数值很低，但不接近于0，通过改进遥感影像暗通道分布从而有效去除薄云。HOT法是利用遥感影像红光波段和蓝光波段相关性高的特点寻找合适的晴空线，用像元到晴空线的距离表征该点云量进而有效去云。He等[11]提出了AHOT(Advanced HOT)算法，先利用空间信息剔除易混淆地物,提高云检测精度，再考虑大气散射乘性效应，提出基于AHOT的虚拟云点(Virtual Cloud Point, VCP)技术扩大由薄云污染降低的数字量化值(Digital Number, DN)方差从而有效去云。Makarau等[12]利用加性薄云污染模型，在传统的暗目标减法(Dark Object Subtraction, DOS)基础上[13]，通过局部搜索整个影像的暗目标构建影像的薄云厚度分布(Haze Thickness Map, HTM)，最后用原降质影像减去薄云厚度分布得到无云影像。这些去云算法大多关注去云后影像细节的锐化、地物间对比度的增强，如何在去云的同时保证数据保真是一项具有挑战性的问题。

近年来，已有学者提出了一些考虑保真的薄云去除算法[14-16]。Shen等[14]基于同态滤波法半自动化地确定不同波段的截止频率，对薄云像素和无云像素分别处理从而获得高保真的无云影像。Singh等[15]基于暗通道先验法提出改进的无云影像恢复算法，通过改进联合三边滤波器以提高大气光的精细估计效果，并重新定义透射率从而降低无云影像的颜色失真程度。Liu等[16]基于加性薄云污染模型，将原始薄云厚度分布先经过纹理去除，再进行大面积地物反射抑制，然后用导向滤波器进行平滑，最后用线性回归求取不同波段的薄云厚度从而得到无云影像。文献[16]算法的去云结果整体具有保真效果，但是依旧存在无云区域过度校正，以及影像中蓝色地物色彩严重失真的问题。

因此，本文在文献[16]算法的基础上提出了一种改进的薄云去除方法，即基于薄云厚度分布(HTM)评估的高保真薄云去除方法。该方法首先利用基于暗通道先验薄云去除算法求得影像的透射率分布图[10]，根据透射率分布图中像素数值接近于1来自动确定无云区域的范围。然后,在HTM中，求取无云区域的平均值,将整个HTM减去无云区域平均值，使得HTM满足无云区域的薄云厚度接近于零，以确保无云区域的保真度。针对蓝色地物区域云量被过估计，造成去云结果中色彩失真的问题，本文利用蓝色地物RGB波段的光谱特性，结合阈值法检测出蓝色地物，用蓝色地物附近的HTM平均值作为蓝色地物区域的云量估计。最后,用降质影像减去不同波段优化的HTM得到最终去云结果。

通过多幅影像的实验发现，与文献[14]算法去云结果、文献[15]算法去云结果，以及文献[16]算法去云结果相比，本文方法的结果影像去云效果良好，同时具有数据高保真的优点。

1 加性模型薄云去除

1.1 云成像模型

薄云污染可以认为是加性污染[12，16],并且薄云污染对不同波段的影响不同[16-17],所以云成像模型可以描述为：

Ic(x)=Jc(x)+Hc(x);c∈{R,G,B}

(1)

其中：x代表像素的位置；I代表降质影像；J代表无云影像；H代表薄云厚度分布；c是{R,G,B}中的任一波段。

1.2 薄云估计

在Makarau等[12]的研究中，薄云厚度分布(HTM)估计如下：

(2)

其中：Ib(x)、Ig(x)分别是降质影像的蓝光波段和绿光波段；W(x)是以x为中心、大小为w×w像素的窗口。本文实验中该窗口大小取3×3像素。

Liu等[16]发现由式(2)得到的Hraw(x)存在将地物反射光作为薄云组成成分，造成降质影像被过度矫正的问题。为了减少地物反射光的影响，文献[16]对Hraw(x)进一步改进：首先为了去除Hraw(x)中的地物纹理，文献[16]采用了可以去除无用纹理和细节的全变差(Total Variation, TV)正则化图像修复模型[18-19]得到平滑的薄云厚度分布Hs(x)；接着,为解决大面积高亮地物的薄云厚度过估计问题，文献[16]用Cannny算子检测出高亮地物边缘，用局部窗口对Hs(x)中高亮地物区域作数值抑制得到Hd(x)；然后,应用导向滤波器[20]去除Hd(x)中的块状效应得到Hg(x)；最后,求取不同波段的薄云厚度分布。

Himpr(x)=Hg(x)-minHg

(3)

Hc(x)=kcHimpr(x)；c∈{R,G,B}

(4)

为了去除气溶胶的影响，式(3)用Hg减去其最小值minHg得到最终的薄云厚度估计Himpr(x)。式(4)中：kc是不同波段的云估计强度，它是在Hg低于平均值的像素集与降质影像I的不同波段中低于平均值的像素集的交集中，以Hg的像素为自变量，以影像的不同波段的像素为因变量，通过线性回归得到；Hc(x)是不同波段的薄云厚度估计。得到Hc(x)后，通过式(1)即可得到无云影像J。

图1是根据文献[16]算法得到的实验结果。其中:图1(a)是大小为300×400像素的GF-2降质影像;图1(b)是最终的云厚度分布估计Himpr(x);图1(c)是去云结果。由图1可以看出，Himpr(x)中尽管云厚度分布范围估计比较准确，但是无云区域的云厚度值并不为0或接近于0，同时无云区域的蓝色和白色物体估计的云厚度分布值也很高，这就造成了去云影像的无云区域存在过度矫正和蓝色地物区域严重失真的问题。

图1 文献[16]算法的薄云估计和去云结果Fig.1 Haze estimation and haze removal results of algorithm in literature [16]

2 薄云去除改进

为解决图1去云结果中反映出来的无云区域过度矫正和蓝色地物失真严重问题，基于加性薄云污染模型，本文在文献[16]得到的最终薄云厚度分布Himpr(x)基础上进行改进，在获得良好去云结果的同时确保影像的高保真。本文方法的流程如图2所示。

图2 本文方法流程Fig. 2 Flow chart of the proposed method

2.1 无云区域保真

首先针对无云区域薄云厚度被过估计，造成无云区域数据失真的问题，本文对Himpr(x)作了如下改进：

Hnew(x)=max(Himpr(x)-K,0)

(5)

其中：

(6)

式中：Dclean是Himpr(x)中的无云区域;M是Dclean中的像素个数;K是无云区域的平均值。理论上，Hnew(x)的无云区域DN值接近或等于0，有云区域DN值大于0，为了避免由于Himpr(x)-K导致负数的出现，这里使用了最大值算子max。

无云区域Dclean的确定影响着最终无云影像无云区域的保真度。为了有效确定无云区域Dclean的范围，本文利用基于暗通道先验薄云去除算法[6-7]中的透射率来估计该区域。其中透射率为：

t(x)∈(0,1]

球形工频电场检测传感器采用球形电容检测感应电压的方法来检测电场,一维球形传感器的结构如图1所示,将一个空心金属球壳切割成两部分,通过绝缘物质将两个半球连接在一起,并在两个半球中间并联一个测量电容,两个半球分别构成了传感器的两个电极[6-7]。

(7)

2.2 蓝色地物云量估计

对于降质影像中的蓝色地物，由于像素在蓝光波段的亮度值远大于绿光波段，因此根据式(1)得到的薄云厚度分布Hraw(x)中容易造成蓝色地物薄云厚度被过估计，最终导致去云结果中蓝色地物色彩失真。

通过大量实验数据发现，在优化后的薄云厚度分布Hnew(x)中，白色地物和蓝色地物区域的DN值都很高。因此，本文先利用阈值选定白色区域和蓝色区域，再根据原始影像中蓝色地物蓝光波段像素值远大于红光波段像素值的特点[21]，剔除白色地物区域最终得到蓝色地物的像素集：

P1={x|Hnew(x)≥T}

(8)

P2={x|Ib(x)-Ir(x)>0}

(9)

Pf=P1∩P2

(10)

其中：x代表像素坐标位置；P1是包含蓝色地物和白色地物的像素集；P2是包含偏蓝地物的像素集；Pf是蓝色地物像素集。通过阈值T可以选定白色地物区域和蓝色地物区域，阈值T的选定不能过大，否则影像中的蓝色地物区域不能全部选中，阈值T也不能过小，否则容易包含非蓝色地物。本文中阈值T取Hnew(x)中最大像素值的70%，即T=max(Hnew(x))*0.7。Ib(x)、Ir(x)分别是原始影像的蓝光波段和红光波段。由于薄云的分布一般是连续和平缓的，因此本文将蓝色地物区域的薄云厚度用周围区域薄云厚度的均值代替，其表达式为：

(11)

其中:Hblue是蓝色地物薄云厚度;Hnew(x)是优化的薄云厚度分布;Dnear是在Hnew(x)中蓝色地物附近人工选取出的薄云区域;n是Dnear中像素点个数。

经过以上的优化，可得到完整的薄云厚度分布Htotal(x)：

(12)

为了去除Htotal(x)中的块状效应，本文应用导向滤波器[20]将其平滑得到最终的薄云厚度分布Hfinal(x)。

本文方法的每一个优化步骤结果如图3所示。其中：图3(b)是根据透射率得到的无云区域二值图像，白色表示无云区域。图3(c)是优化后的薄云厚度分布Hnew(x)。可以看出，Hnew(x)在无云区域的DN值接近或等于0，这能保证恢复的无云影像在无云区域的高保真度。图3(d)是蓝色地物的检测结果, 可以发现，图中的蓝色地物位置检测比较准确。图3(e)是最终的薄云厚度分布Hfinal(x)。与图3(c)相比，图3(e)的Hfinal(x)在蓝色地物区域的DN值相比Hnew(x)小很多，可以避免蓝色地物云量被过估计。

图3 本文方法各步骤结果Fig. 3 Results of each step of the proposed method

2.3 高保真无云影像

通过逐步优化薄云分布，单独估计蓝色地物云量，最终得到的无云影像J表示如下：

Jc(x)=Ic(x)-kcHfinal(x)

(13)

其中：I是降质影像；kc是根据文献[16]算法中通过线性回归得到的不同波段云估计强度；Hfinal(x)是最终优化的薄云厚度分布。图3(f)是本文改进方法得到的无云影像，对比图3(a)的降质影像可以发现，图3(f)的无云影像不但具有良好的去云效果，而且无云区域的地物色彩基本不变，蓝色地物色彩自然，即本文方法得到的无云云影像具有高保真的特点。

3 实验结果及分析

本文采用三幅大小均是512×512遥感影像作为实验对象。图4是分辨率为30 m的Landsat7遥感影像；图5是分辨率为5.8 m的ZY3遥感影像；图6是分辨率为2.0 m的WV-2遥感影像。数据包括了城市、农田、草地等多种地物类型。本文采用文献[14]方法、文献[15]方法以及文献[16]方法与本文方法进行比较。

图4 不同方法去云结果比较(Landsat7遥感影像)Fig. 4 Haze removal result comparison of different methods (Landsat7 remote sensing image)

图5 不同方法去云结果比较(ZY3遥感影像)Fig. 5 Haze removal result comparison of different methods (ZY3 remote sensing image)

从图4～6可看出，4种算法都具有良好的去云效果。文献[14]去云方法会使降质影像的云区变得有些模糊，使得地物细节不够清晰，且去云影像的无云区边界色彩过渡不自然。文献[15]去云方法对影像地物数据的保真度较差，偏绿地物太过鲜艳，色彩不自然。文献[16]去云方法对非蓝色地物整体有一定的保真效果，然而在蓝色地物区域(图5(d)和图6(d)中方块)去云后色彩与原图中的蓝色地物(图5(a)和图6(a)中方块)相比过于鲜艳，失真严重。本文方法得到的去云结果不但在无云区域有较高的保真度，而且影像中的蓝色地物(图5(e)和图6(e)中方块)去云后色彩自然，达到高保真去云目的。

为了定量评估本文方法的去云效果，采用广义奇异值分解(Generalized Singular Value Decomposition, GSVD)作为评价因子[22]。GSVD是用于评估图像受到模糊或噪声影响的程度,薄云覆盖是光学遥感影像被模糊的主要原因。GSVD值越小，则影像的去云效果越好。

不同方法的GSVD统计值如表1所示。从表1可看出，四种去云方法中，本文方法对应的三波段GSVD平均值总体较小，表明本文方法具有理想的去云效果。

图6 不同方法去云结果比较(WV-2遥感影像)Fig. 6 Haze removal result comparison of different methods (WV-2 remote sensing image)表1 不同方法的GSVD统计值对比Tab. 1 Comparison of GSVD statistical values of different methods

图像波段降质影像文献[14]方法文献[15]方法文献[16]方法本文方法图4蓝波段25.8524.6225.1824.7925.14绿波段24.5824.2323.9824.0424.08红波段23.3523.5623.3123.2123.20平均值24.6024.1424.1624.0124.14图5蓝波段22.5120.8118.6719.1812.15绿波段17.3514.8615.4716.9317.30红波段19.6819.1518.5019.6919.44平均值19.8518.2717.5518.6016.30图6蓝波段15.1713.6313.7613.9514.35绿波段13.2913.7212.8113.0912.37红波段10.029.459.709.9810.01平均值12.8312.2712.0912.3412.24

为了定量评价降质影像无云区域的保真度，本文采用平均相对误差(Mean Relative Error, MRE)作为评价指标，结果记作δ:

(14)

其中：Dclea是降质影像I(x)和恢复的无云影像J(x)的同一片无云区域；m是无云区内像素点个数。δ值越小保真效果越好，该评价指标MRE的统计值如表2所示。

表2 不同方法MRE统计值对比Tab. 2 Comparison of MRE statistical values of different methods

由表2可看出，四种去云方法都具有高度的保真特性，文献[15]方法去云结果的数据保真性能相比其他方法略差一些。本文方法的无云区域保真性能最好。综合两个评价指标可以看出，本文方法能得到高保真的最优去云结果。

通过对不同数据范围、不同分辨率大小、多种地物类型的薄云污染影像数据实验得出，本文方法得到的无云影像都具有数据高保真效果，表明本文方法对不同薄云污染影像去云具有高度适用性。

4 结语

本文提出了一种基于薄云厚度评估的高保真去薄云方法，不仅能够满足无云影像中无云区域数据的高保真度，还能够解决传统算法蓝色地物色彩失真问题。针对影响算法的两大因素：无云区域Dclean的确定和阈值T的选取，本文通过透射率分布图中像素数值接近于1实现了自动化确定无云区域的范围，然而在确定蓝色区域的过程中，本文是以经验阈值T选取蓝色地物和白色地物区域，自适应求取最优阈值T将是我们后续的改进工作。