改进暗通道先验的高分二号影像快速去薄云

潘勇卓，谢洪斌，杨雪，席亚文，戚伟迅

(1.重庆地质矿产研究院 外生成矿与矿山环境重庆市重点实验室，重庆 400042；2.煤炭资源与安全开采国家重点实验室重庆研究中心，重庆 400042)

0 引言

近年来高分系列卫星在国家高分专项的推进下，数据质量不断提升，尤其是高分二号影像空间分辨率可达0.8 m，突破了亚米级的重要门槛，目前已被国防、科研、建设等各行各业广泛使用。然而高分二号卫星回访周期相对较长，数据量十分有限，所获取的信息质量依赖于地表的光谱反射。通常云层也会作为一种地表物质记录在遥感影像上，导致原本地表信息被遮挡，降低了影像的利用率和可判读性。随着高分系列卫星影像分辨率的不断提高，云层影响变得更加不容忽视[1]。云层可大致分为厚云和薄云，厚云遮挡了绝大部分的地面辐射，采用其他时相的影像进行插补是唯一的解决方法[2]；而薄云仅削弱了地物的部分反射特性，因此可以基于单幅影像的特征降低薄云干扰，重建下垫面地物信息，提高高分影像的利用率[3]。

目前，国内外单幅影像去薄云的方法主要分为频域去云和空域去云两类。前者专门针对图像上薄云雾进行滤波处理，后者根据多光谱图像不同波段对云层的反射率差异进行云层检测和处理[4]。频域去云领域较为成熟的方法有同态滤波[5]、小波变换[6]、自适应滤波[7]等。频域去云基本可达到薄云去除的效果，但通常会损失部分地物的低频信息，导致模糊和纹理缺失。关于空域去云算法的研究较多，其中，直方图匹配法[8]是对云区影像和晴空区影像的地类进行直方图匹配，但去云结果图可能出现边缘过渡带；最大化局部对比度法[9]通过调整图像的对比度来实现薄云去除，但对饱和度较高的原图像易造成色彩失真的效果；基于HOT[10]的去薄云方法利用红蓝波段的高度相关性生成HOT图，根据HOT图反映的云雾影响程度分层实施暗目标减法，该方法在植被覆盖区取得了很好的去薄云效果，但对水体、裸地及人工地物等地物类型十分敏感，容易出现过度校正的情形；暗通道先验法[11]是近年来计算机视图领域提出的去雾算法，该方法具有适用性广和速度快的优势。随后一些研究者基于该方法进行了一些改进处理，如通过设置阈值来提高透射率图的计算速度[12]、通过导向滤波优化透射率图的精细程度[13]、通过设置暗通道权值完善图像不连续处的处理效果[14]、将RBG空间转换到HIS空间并对图像进行分割从而提高图像的清晰度[15]、对景深边缘和非景深边缘采用不同模板进行滤波从而降低边缘的不连续效应[16]等。这些改进的暗通道先验法主要针对RGB普通光学照片去雾，遥感影像上的“薄云”与这些“雾”类似，近年来相关学者将该方法应用于遥感影像，获得了较好的去云雾效果[17]。但该方法在处理成片的高亮区域时将出现过度校正的情形，通常在遥感影像上体现为厚云过校。此外，由于蓝光波段对云层最为敏感，直接应用暗通道先验法对遥感影像进行处理将造成色彩失真的问题。

本文对暗通道先验法[11]进行了改进，提出了一种适用于高分二号影像的薄云快速去除方法。该方法通过设置最大大气光阈值，有效改善厚云及其他高亮地物的过校问题；计算透射率分布图时，对原图进行下采样处理，以提高透射率的计算速度；利用暗通道法计算去云图像时，非独立地计算各波段的改正量，而是结合高分二号影像各波段间散射关系，以某一波段为起始波段，通过波段间散射强度的校正系数确定其他波段的改正量，从而使去薄云影像的各波段间关系更符合大气散射模型，避免了RGB合成色彩失真的问题。

1 实验数据说明

本文使用的高分二号影像获取时间为2016年6—8月，覆盖区域为重庆市某区的局部范围，空间分辨率为0.8 m，包含蓝、绿、红和近红外波段。为便于处理，本文将4个波段合成为.tif格式文件作为实验输入。

2 原始暗通道先验法

2.1 暗通道先验理论

暗通道先验理论由He等[11]提出，是一种基于对大量自然无雾图像的观察统计而得到的计算机视图规律：在绝大多数非天空、无雾彩色图像的局部区域里，某些像素总会有至少一个颜色通道具有非常低甚至为零的灰度值，这样的像素成为暗像元。暗像元主要存在于以下情况中：①植被、建筑、山体、车辆等各种地物的阴影；②绿色、蓝色、红色等某一颜色通道亮度值极高，而某一颜色通道亮度值较低的彩色物体；③黑色的树干、人工地物、裸地等。暗通道图像可通过式(1)来描述：

(1)

式中：c表示r、g、b中的一个通道；Jc是图像J的3个通道中某一个通道的灰度值；Ω(x)是以像元x为中心的局部领域区域；Jdark为暗像元的灰度值，由暗通道先验理论可知，对于图像的无云雾区域Jdark→0，而云雾区域的暗像元灰度值会变高。

2.2 云雾图像退化模型

在计算机视觉领域，Koschemieder模型方程[18]所描述的云雾图像退化模型被广泛使用：

I(x)=J(x)t(x)+A(1-t(x))

(2)

式中：I(x)为观测图像的辐射强度，代表待去除云雾的原始图像；J(x)为来自目标场景的辐射信息，代表去除云雾后的退化图像；t(x)为场景辐射在大气传输过程中的透射率，代表透射率分布图；A为大气光值。该式第一项J(x)t(x)描述了目标场景辐射信息在介质中的直接衰减，第二项A(1-t(x))代表大气光经过衰减后导致的场景颜色偏移。文献[11]指出，当大气条件均匀且同质时，透射率t(x)可表示为：

t(x)=e-βd(x)

(3)

式中：β为大气散射系数；d(x)为场景深度。该式表明场景辐射的透射率与云雾浓度和景深成指数衰减。与普通室外照片相比，高分辨率遥感影像中云雾的分布更为均匀，且地表到卫星传感器之间的距离相对固定，因此遥感影像的透射率分布图t(x)更为平滑。

图像去云雾的目的，就是求取A和t，从而通过原始图像I还原得到退化图像J。若图像I共包含N个像素，则云雾图像的退化过程将包含RGB三通道的3N个约束方程和4N+3个未知数，这显然是一个无解的欠约束问题，需要增加约束条件，而暗通道先验理论就为其提供了一个非常好的约束条件。

2.3 暗通道先验去云雾流程

1)估算大气光值。由式(1)可知，暗通道图像反应了原始图像中的云雾分布情况，暗通道值越高，云雾越浓厚。单幅带云雾的图像中，大气光值A可近似为最浓云雾区域的灰度值，因此借助暗通道图即可从原始图像中估算该值。He等[11]设计了A的估算方法：从暗通道图中按照灰度的大小取前0.1%的像素，在这些像素对应的位置中，从原始图像取灰度的最大值作为大气光值A的估计值。

(4)

根据暗通道先验理论[11]，由式(1)可推导出：

(5)

将式(5)带入式(4)可得到：

(6)

少量云雾的存在可为图像保留一定的透视深度感，使其更真实，因此引入控制参数ω(0<ω≤1)来保留少量云雾。透视率的粗略估计式为：

(7)

3)优化透射率分布图。实际影像中，云雾对物体的影响程度分布是连续的、平滑的，而式(7)估算的透射率是在局部领域区域Ω(x)内大气透射率保持一致的假设下计算得到的，所形成的透射率图将出现斑块效应，与实际不符。为了获得更精细的透射率图，He等[11]首先提出了Soft Matting方法进行内插优化，消除了透射率图的斑块效应，但该方法的最大缺点在于处理速度较慢。随后He等提出了导向滤波[19]的方式来改进透射率图的处理，该方法主要集中于简单的方框模糊，不仅能大幅提高处理速度，而且使透射率图的效果更精细。导向滤波指出，输出图像的轮廓由输入图像决定，而输出图像对场景结构信息的保持性由指导图像决定。导向滤波的模型为：

Qi=akIi-bk，∀i∈wk

(8)

式中：Qi为输出图像；Ii为指导图像；ak和bk为窗口wk中与输入图像有关的线性常量系数。对该式求导变形为：

(9)

该模型保持了指导图像和输出图像边缘变化的一致性，而输入图像决定了输出图像的轮廓，因此可通过最小化输入图像与输出图像间的差异来计算线性常量系数ak和bk。定义能量常数为：

(10)

4)还原退化图像。在获得大气光值A和透射率分布图t后，通过式(2)变形，即可还原云雾退化图像。而在实际应用中，当透射率t很小时，会导致退化图像整体向白场过度。He等[11]通过设置最小阈值t0对透射率进行控制，最终得到的图像还原计算式为：

(11)

3 改进的暗通道先验法

3.1 最大大气光阈值法

暗通道先验法直接用于遥感影像的处理时，将会出现厚云和其他高亮地物的过度校正问题。这是由于该方法中大气光值A最终选取的是原始图像中的某一个像素点的值，则各通道的A值很有可能全部接近255，从而导致处理后的图像偏色和出现大量色斑。

3.2 基于下采样方法的透射率快速计算

He等[11]利用导向滤波对透射率的计算过程进行了优化，使透射率图的处理精度得到了有效提高。但高分二号影像分辨率可达0.8 m，单景影像所包含的浮点计算量极大，直接使用原始图像计算将导致透射率图的处理时间太长。考虑到云雾的分布是连续、成片的，对于遥感影像薄云去除，并不需要米级的参数矫正，因此可以通过合并像素降低分辨率的方式近似估算原图的透射率图，以达到降低计算量而又不影响校正结果的目的。为此，本文采用下采样的方法提高透射率的计算时效，具体实现方式为：

①在求取透射率之前，首先按照采样率r对原图进行下采样操作，将其尺寸缩小；

②按照原算法，求取下采样后图像的透射率；

③通过插值的方法上采样得到原图的透射率。

下采样操作略微降低了透射率图的细腻程度，而计算量得以有效减少。为探讨不同采样率r对暗通道先验法去薄云结果的影响，本文从高分二号卫星影像中截取局部含云区域(截图大小1 024×1 024)进行实验。图1展示了r分别取0.25、0.5、0.75和1(即未采用下采样法)时的实验结果，图2统计了r在[0.2，1]范围内去薄云图像均值和方差的变化情况，图3展示了r对处理速度的影响。

图1 不同采样率下薄云去除的效果

图2 不同采样率下图像均值和方差的变化

由图1可以看出，不同采样率的暗通道先验法去薄云效果差异不大，但随着r降低，图像略微变暗。由图2可见下采样方法导致图像均值逐渐降低，而方差变化并不大，说明该方法的使用对去薄云的整体成像效果影响较小；而随着r的进一步降低，图像均值的降低程度增大，当r<0.25时，经过下采样操作后的图像均值相比原始方法降低了15%以上。因此，为合理控制图像的视觉效果，可选择采样率的取值区间为[0.25，1)。

根据图3，未采用下采样操作时，原始暗通道先验法耗时为51.31 ms，采用下采样操作后，随着r的降低，图像处理耗时随之减少，说明下采样方法有效降低了透射率图的计算量，因而去薄云效率得以提高。当r<0.7时，下采样操作使得处理速度提升了1倍以上；当r逐渐增大至0.9时，经过下、上两次采样后处理速度的提升将会逐渐减小。值得注意的是，当r取值更大时，处理耗时将比原始方法更多，如r=0.99时耗时达59.82 ms，此时下采样操作不具有收益性。因此，为有效提升处理速度，可选择采样率的取值区间为[0.2，0.7]。

图3 不同采样率下处理速度的变化

综上，合理的下采样操作将明显提升处理速度，且成图的视觉效果与原始算法无太大差别。当采样率在[0.25，0.7]区间范围内取值时，下采样方法将具有较明显的收益性。

3.3 考虑波段间散射关系的云雾退化过程

不同于普通光学照片，遥感影像中各波段对云层的反应存在差异：遥感影像中，受到云层影响时，各波段的DN值都会升高，但蓝色波段对云的干扰更敏感，云越浓厚，DN值升高的幅度越大[3]。不同像素的云雾覆盖情况对波段DN值的影响存在差异，而暗通道先验法应用于遥感影像时，相互独立地寻找每个波段的还原值，致使像素内不同波段的值偏离于真实的大气散射模型，甚至导致影像某一(或几个)波段被过度校正，出现RGB合成影像色彩失真的问题。

为此，本文并不独立地寻找各个波段的还原值，而是考虑大气散射模型，以某一波段为起始波段，根据各波段间散射强度的比值关系设置校正系数，推算其余波段的散射改正量，从而约束波段被过度矫正的问题。根据瑞利散射规律，散射光强度ADN与入射光波长λ成反比关系，而在薄云条件下该反比关系可限制在λ-0.7至λ-0.5之间[20]。

去薄云过程中，首先选择对云层最敏感的蓝光波段为起始波段，通过暗通道先验法求取还原值，并以蓝光原始值DNB与还原值DN′B之差作为起始改正量，而其余波段的散射改正量将由波段间的校正系数直接产生，进而求取其余波段的还原值。本文选取λ-0.7关系为波段间校正系数的计算参考，高分二号卫星各波段参数及以蓝光为起始波段产生的校正系数如表1所示。

表1 高分二号卫星参数及波段校正系数

首先通过暗通道先验法求取蓝光波段还原值，则蓝光波段的散射改正量为：

ADNB=DNB-DN′B

(12)

根据表1的校正系数，以蓝波段为起始波段计算绿、红、近红外波段的散射改正量ADNG、ADNR、ADNNIR分别为：

(13)

由此可求得绿、红、近红外波段的最终还原值为：

(14)

4 高分二号影像薄云去除

4.1 改进的暗通道先验去薄云流程

根据改进的暗通道先验法，高分二号影像的薄云去除流程如图4所示。经实验，本文以目视效果最佳为原则相关设置以下参数：暗通道变换最小窗口大小为15×15，云雾控制参数ω=0.95，透射率的最小阈值t0=0.1；以原图的灰度图作为导向滤波的导向图，滤波窗口大小为60×60；最大大气光阈值A0=220；下采样率为0.25。其余参数设置已在前文中给出。

图4 改进的暗通道先验法去薄云流程图

4.2 实验结果与对比分析

根据上述改进的暗通道先验法，截取含云的高分二号影像进行去薄云实验，并分别与空域去云领域的HOT算法、原始暗通道先验法的去薄云影像进行对比，同时采用无云或少云时相的高分一号卫星影像作为客观参考，从校正后影像的整体直观效果对比(图5)、局部纹理色彩对比(图6)和云雾校正真实性统计(表2)等方面综合评价各方法的去薄云效果，最后评价处理效率(表3)。其中云雾校正真实性统计是通过将校正后影像的含云区域与参考影像对比、无云区域与原始影像对比，定量统计各方法形成的图像失真面积来实现的。

图5 去薄云效果对比图

图6 局部细节对比图

1)去薄云效果评价。从直观对比可以看出，HOT法去薄云效果似乎最为明显(图5(b)、图5(g))，但整体色彩偏绿。通过局部放大可发现，该方法降低了影像细节的真实性，且存在严重的色彩失真现象。其中，图6(b)中厚云中心区域被过度校正，形成了黑色斑块；图6(g)中稍厚的云层被“绿化”为与植被相似的地物，使得整体图像存在大量“云层下的植被被还原”的误导信息，校正后显示为虚假的草地，这与该

表2 云雾校正真实性统计

表3 各类方法处理时间表

区域无云影像中(图6(j))原有的农田与房屋信息严重不符。同时大量蓝色建筑同样被“绿化”(图6(l)、图6(q))，原有房屋难以从色彩上判别，甚至某些高亮的蓝色建筑色彩信息直接丢失。同时裸地、硬化地面等红棕色、灰色地物也出现泛白的情况，去薄云后影像色彩效果较差。通过表2可以看出，该方法导致的失真面积达1.258 km2，占两幅影像总面积的3.558%，意味着校正后影像将丢失这部分真实信息。此外，在局部图中还可以发现，该方法导致大量噪声，纹理效果大大降低，不利于影像后期使用。综上，HOT法去薄云的纹理色彩、校正真实性均较差，整体校正效果与参考影像差距较大。

原始的暗通道先验法对于薄云有较明显的去除效果(图5(c)、图5(h))，而稍厚的云层则予以保留，且在图像纹理的处理上更贴近原始图像(图6(m)、图6(r))，避免了色块的噪声影响。同时该方法的去薄云图像比HOT法更符合事实，薄云区域地物还原情况更接近参考影像。但该方法使整体图像偏蓝、偏暗，尤其是稍厚的云层“蓝化”现象严重(图6(h))，且在亮度较高的厚云区域仍存在过度校正的问题(图6(c))。经统计，该方法云雾校正的失真面积约为0.109 km2，占2幅影像总面积的0.308%，去薄云后影像的真实性高于HOT法，但色彩效果仍与参考影像存在一定差距。

本文方法在去薄云的效果上继承了原始暗通道先验法的优点，而对其存在的问题进行了明显的优化。从整体效果上看，本文方法去薄云后的图像(图5(d)、图5(i))仍保留了部分厚云，色彩更接近原始影像和参考影像，克服了原始暗通道先验法的偏蓝问题；从局部细节上看，本文方法避免了稍厚云层的蓝化现象和高亮厚云的过校问题(图6(d)、图6(i))，效果图中未发现严重的失真现象，在3种算法中保持了最好的云雾校正真实性。此外，该方法对图像纹理、色彩的还原都更为逼真(图6(n)、图6(s))，细节上与参考影像中的无云情况最为接近。

2)处理速度评价。为展示本文方法的速度优化，表3给出了本文方法与原始暗通道先验法在图像处理时间上的对比(时间结果为3次计算取平均)。由于HOT法为非自动方法，因此不参与处理速度评价。本次实验采用Windows10系统，算法仿真采用Visual Studio2013平台，计算机处理器为Core i7-4790 CPU@3.60 GHz，内存8 GB。通过表2可以看出，本文方法的去薄云处理时间约仅为原始暗通道先验法的1/6，处理速度明显更优。

综上实验结果表明，本文方法能较好地降低薄云对高分二号影像的干扰，成果图像的纹理色彩、校正真实性和处理速度相比HOT方法和原始暗通道先验法都最具优势。

5 结束语

本文将暗通道先验法引入到遥感影像的处理上展开研究，并结合卫星波段特点，提出了一种适用于高分二号影像薄云快速去除的改进方法。改进之处总结为以下3个方面：

①针对暗原色先验法对厚云等高亮区域的过校问题，提出了最大大气光阈值的解决策略，避免了该类区域的失真；

②在原有透射率优化方法的基础上，提出了下采样和插值算法简化输入对象，大幅提高了处理速度；

③针对成果图像的色彩偏蓝现象，结合高分二号传感器波段间散射强度相互关系，以蓝光为起始波段推算其余波段，保证不同波段的改正量符合相对散射模型，从而使成果图像的色彩更贴合实际。

综上所述，本文方法兼顾了去薄云图像的纹理色彩、校正真实性和处理效率，对于单幅高分二号影像薄云快速去除具有较高的实用性。当然，本文方法尚有待改进之处需要进一步研究，如对其他卫星影像适用性的探讨、如何建立可靠性更高的评价体系来衡量去薄云后影像的质量等。