基于边缘敏感递归滤波的彩色航拍图像云检测

廖 斌，付忠旺

（湖北大学计算机与信息工程学院，湖北武汉430062）

基于边缘敏感递归滤波的彩色航拍图像云检测

廖 斌，付忠旺

（湖北大学计算机与信息工程学院，湖北武汉430062）

为了有效检测出彩色航拍图像中的云区，提出一种基于边缘敏感递归滤波的彩色航拍图像云检测方法。通过基于自适应引导权的边缘敏感递归滤波算子构造航拍图像的差异映射，并结合基于类内方差的递归阈值选择方法，将输入的航拍图像区分为待定云区和非云区。利用边缘敏感递归滤波算子对原始航拍图像累进滤波，有效抽取分层的图像细节，构造航拍图像的细节映射以精化待定云区。利用客观性能指标对检测结果进行评价。实验结果表明，相较于其他相关研究工作，提出的方法有效易行，计算复杂度低。

图像处理；图像分类；云检测；彩色航拍图像；边缘敏感递归滤波

0 引 言

作为一类典型的空间信息数据源，航拍图像以其成像清晰度高，信息客观丰富以及时效性强等优点而在天气分析、灾害与环境监测、资源勘探等诸多领域有着广泛应用。但是，由于云覆盖了大约一半的地球表面，航拍图像中通常含有云区，极大影响了图像识别、图像分类等处理的精度，使得数据的利用率大幅度降低。航拍图像的云检测对于缓解海量航拍图像数据的传输、存储、分类与检索的压力，提高数据利用率具有非常重要的意义。如果仅仅通过人工交互完成云检测，会非常耗时且检测精度也难以保证。因此，有效自动的云检测成为航拍图像处理与应用的首要任务及关键技术之一。______

已有的云检测方法大多利用多光谱卫星图像融合技术，根据云在不同谱段下的反照率，选取合适的阈值进行云检测［1］。随着MODIS等高分辨率探测器的出现［2］，国内外学者运用多光谱高分辨率图像提高了云区与非云区的判别精确度。不同于多光谱卫星图像，彩色航拍图像只含有红色-绿色-蓝色（red-green-blue，RGB）三通道的光学数据，没有其他任何辅助通道信息。由于基于多光谱数据的云检测方法大都高度依赖RGB波段之外的可用频谱波段，因此不能直接应用于彩色航拍图像的云检测。文献［3］基于双边滤波算子提出一种彩色航拍图像的云检测框架。但是，双边滤波算子不适用于对任意尺度的细节提取，当增大平滑参数多次迭代时，其不能很好地保持边缘，会导致梯度反转问题［4］。同时，双边滤波算子对于边界模糊的情况也不能较好地处理。因此，该框架会导致检测出来的非云区含有误认为云区的细微聚类，且检测出来的云区边界不理想，云区内部还会出现洞隙。这些问题都需要借助于后续处理来解决。此外，该框架还要进行多尺度双边滤波计算，非常耗时。尽管已有一些加速双边滤波计算的方法［56］，但是这些方法通常采用量化或者下采样的策略，不可避免的对计算结果产生影响。

彩色航拍图像的云检测本质上是一个图像分割的问题。虽然，图像分割领域已经有了大量的研究工作，但并没有形成通用的分割理论。现有的分割算法多是解决某类具体问题。自动图像分割方法［79］通常无法有效处理复杂的彩色航拍图像。而交互式图像分割方法［1011］都需要大量的人工交互，无法直接处理大数据量的彩色航拍图像。

边缘敏感的图像处理方法在图像与视频编辑、颜色表征感知等领域都有广泛的应用。最流行的边缘敏感滤波方法是双边滤波，但是其执行效率较低，对于模糊边缘无法有效感知。与之相比，引导滤波在处理图像边缘时表现的更好［12］。然而，其用到了局部线性模型，当滤波的核半径过大时，滤波结果仍然会出现边缘模糊的问题。同时，引导滤波对整幅图像都采用统一的滤波半径，对于边缘、细节等特征难以有效的感知。文献［13］提出基于域变换的线性递归滤波算法，大大降低了滤波的复杂度，可以实时地对图像和视频进行保边缘滤波。但是，该算子仅利用像素间颜色特征的差值来定义颜色特征权，在处理边界模糊的复杂边缘时，仍然难以有效地保边缘，同样也会出现梯度反转的问题。

基于以上分析，本文提出一种有效的彩色航拍图像边缘敏感云检测方法。研制一种基于自适应引导权的边缘敏感递归滤波算子。该算子不仅考虑了航拍图像的空间特征信息，同时也考虑了颜色特征信息。利用边缘敏感递归滤波算子构造航拍图像的差异映射，有效地拉伸云区和非云区的差异。通过基于类内方差的递归阈值选择方法，结合差异映射，将输入的航拍图像区分为待定云区和非云区。利用边缘敏感递归滤波算子对原始航拍图像进行累进滤波，有效获取分层的图像细节，并构造细节映射，用以精化待定云区，去除其中的非云区，从而得到最终的云检测结果。由于边缘敏感递归滤波算子天然适于进行图形处理器（graphics processing unit，GPU）加速计算［14］，本文方法从根本上相比其他方法具有更好的GPU并行执行效果，能够在获取满意的云检测结果的同时保有更高的性能。

1 边缘敏感云检测

1.1 基于自适应引导权的边缘敏感递归滤波

设Xn为输入的离散信号，Jn为输出的离散信号，n是离散信号的索引，传统的一维递归滤波算子可定义如下：

式中，a∈［0，1］是反馈系数［13］，可定义如下：式中，σH为空间滤波参数，本文采用文献［13］中的推荐取值σH＝15。

由于该算子是一维线性的，只能沿着输入信号的每一维度执行。并且，该算子的输出依赖于输入和先前的输出信号，所以其是非对称的。因此，其必须沿着输入信号的每一维度双向执行，才能完成一次滤波处理。图像可以视为二维离散信号。利用该算子对图像滤波时，应该沿着图像的水平方向从左向右执行，再从右向左执行。然后，再沿着竖直方向双向执行。

同时，由于a只是在信号的空间域上定义，即该算子没有考虑信号的范围信息，导致其会产生无限脉冲响应。如果直接采用该算子对彩色航拍图像进行滤波处理，将无法保持图像的边缘信息，会引起边缘的模糊。受文献［12- 13］的启发，本文拟将引导权引入到该算子。但是，传统的引导权采用统一的滤波半径，难以精确地感知彩色航拍图像边缘和细节特征。基于此，提出一种基于自适应引导权的边缘敏感递归滤波算子，定义如下：

式中，Wn，n－1为自适应引导权，定义如下：

在式（3）中，Y表示引导图像，其可以与原始图像一致，也可以是不同的图像。本文设定Y与X一致。式（4）中，ωn、ωk表示在Y内中心像素位于索引n、k且半径为rn、rk的滤波窗口。其中，如果rn、rk取值过大，滤波结果会出现边缘模糊的问题。如果rn、rk取值过小，则会导致滤波窗口中的邻域信息不够充分。本文以文献［12］中的窗口半径推荐取值60为基准，经过实验测算且不失一般性，取rn＝80，rk＝60。｜ωn｜、｜ωk｜分别为ωn、ωk中所包含像素点的个数。μk和分别为ωk中所包含像素点均值和方差。ε是规则化参数，其主要为了增强计算鲁棒性，避免可能为0导致除式无意义的问题。根据ε≪的原则，本文取ε＝10－6。εn是自适应补偿量参数，受文献［15］的启发，εn定义如下：

式中

由于边缘敏感递归滤波自身具有可并行处理的特性，因此，本文采用了基于统一计算设备架构（compute unified device architecture，CUDA）的GPU计算环境来实现程序的高性能运行。

1.2 待定云区检测

由于色相-饱和度-亮度（hue-saturation-intensity，HSI）模型可以在最大程度上符合人类的视觉感知，因此，首先将输入的彩色航拍图像的RGB模型转换成HSI模型。传统RGB模型转换成HSI模型的变换公式见文献［16］。本文根据大量的观察和统计结果［3］，改进变换公式如下。设

则H、S、I的计算式为

式中，R、G、B、H、S和I分别表示红色、绿色、蓝色、色相、饱和度和亮度。

由于彩色航拍图像中的云区通常具有较高的I和较低的H，受文献［3］的启发，约束I和H的范围为［0，255］，构造差异映射来凸显云区和非云区的差异。定义输入信号为＋255）／（＋255），根据式（3）进行滤波，则输出信号即为差异映射。其中和分别表示X中像素点n的I和H。

将差异映射同样约束到范围［0，255］，则由图1给出其可视化效果，且与文献［3］方法的显著映射进行了对比。可以看到，在图1（c）中，本文的差异映射成功地拉伸了云区和非云区之间的差别，同时很好地保持了边缘。而在图1（b）中，文献［3］方法的显著映射虽然也可以拉伸云区和非云区之间的差别，但是在很多稀薄的云区处未能保边缘，甚至还会丢失云区信息。

基于差异映射，通过选择合适的阈值即可区分出输入图像中的初步待定云区和非云区。通常可以采用Otsu方法［17］来实现这一过程。该方法假设输入图像含有两类像素，计算最优的阀值来区分这两类像素，并且要求两类像素的类间方差取最大值。但是，该方法并没有考虑两类像素的类内方差因素的影响，对于阈值的选择会产生一定的偏差。本文将类内方差作为阈值计算的约束条件，使得到的阈值更合理。同时，鉴于航拍图像往往较为复杂，不宜采用全局阈值选择方案。基于上述分析，提出一种基于类内方差的递归阈值选择方法。选择的阈值T计算如下：

式中

pi为灰度级i的出现概率；λ0ω0σ0＋λ1ω1σ1为类内方差，λ0、λ1为平衡因子。

图1 差异映射与显著映射对比

利用T可将差异映射区分为初步待定云区和非云区。将得到初步待定云区作为一幅新的图像再次通过式（11）计算新的T，再次将差异映射区分为初步待定云区和非云区。如此，递归执行下去，直至满足结束条件时停止，从而获取最终的初步待定云区。

本文将初步待定云区与非云区的类间方差和面积差作为递归运算的结束条件。设第i次阈值选择后，非云区和初步待定云区的面积分别为Ai0、Ai1，则归一化的非云区与初步待定云区的面积差的绝对值为

式中，Ai0＋Ai1是当前整幅图像的面积。此时非云区与初步待定云区的类间方差为

可以看出，经过多次递归迭代，ΔAi将逐渐增大，而σi则将逐渐减小。由此，结束条件参数Gi可定义为

设结束条件阈值为G，如果满足Gi≥G，则结束递归运算。由于G控制递归运算的迭代次数，如果G设置过大，则会将部分云区误认为非云区，导致检测结果的质量下降。如果G设置过小，则部分非云区还未从云区中去除，同样影响检测结果。本文对文献［18］数据库中的数据进行大量测算，选取G＝0.52。

在具体实现中，根据上述过程计算得到的T，结合大量的统计结果［3］，可定义初步待定云区如下：如果差异映射中的像素点n满足当T＜95时，J≥95，或当100≤T≤145时，J≥T，或当T＞145时，J≥145，则满足条件的像素点集合就构成了初步待定云区。

通过上述方法将差异映射区分为初步待定云区和非云区。同时，由于航拍图像中云的H取值范围为［50，80］［3］，需将初步待定云区中不满足条件50≤X≤80的像素点去除，从而得到最终的待定云区U。为了更加直观地呈现U和非云区，可视化定义如下

图2比较了本文方法得到的待定云区和文献［3］方法得到的待定云区的结果对比。首先，从图2（b）和（c）中，可以看到，本文的差异映射相较于文献［3］方法得到的显著映射，获取的云区边缘更加准确，也没有丢失云信息。而在图2（e）中，利用本文方法获得的待定云区，更加符合航拍图像的真实状况。在图2（d），文献［3］方法得到的待定云区则有一些错误检测。

虽然本文方法得到的待定云区非常接近真实的云区，但是，如果输入的航拍图像中含有高亮度的地面对象，本文的递归阈值选择方法也无法从差异映射中将其区分出来，其会被检测为云区。在图2（e）中，可以看到有一些地面对象被误选为云区。因此，需要进一步去除这些非云区，以获取更加理想的检测结果。

图2 待定云区检测结果对比

1.3 待定云区精化

根据上述分析，待定云区中可能会含有一些非云区。这些非云区主要是航拍图像中高亮的地面细节。而真实的云区往往是以相似的像素成团出现，几乎不含图像细节。基于此，利用边缘敏感递归滤波算子抽取这些地面细节，构造细节映射，用以精化待定云区。

首先，利用边缘敏感递归滤波算子对输入的彩色航拍图像Xn进行滤波处理，可以将其分解为基层和细节层。其中，滤波结果Jn即为基层。细节层Dn可定义如下：

如果利用边缘敏感递归滤波算子对Xn进行累进滤波可以得到M个滤波结果｛。根据式（3），迭代得到的过程如下：

式中

为了有效防止噪声以及云的模糊边缘的影响，本文并不基于某一细节层来构造细节映射，而是抽取M个连续的细节层｛，其中

参数M决定着累进滤波的迭代次数。M取值越大，则获取的细节信息越丰富，但同时会导致计算量过大。本文采用文献［4］的推荐取值M＝5，构建5个细节层。由于输入的航拍图像可能会含有较大的噪声，通常不考虑第一个细节层D。然后，对抽取的细节层进行加权平均，从而构造出高质量的细节映射XDetailn：

式中

在图3中，将本文的细节映射与文献［3］方法得到的细节映射进行比较。可以看出，本文的细节映射在保持云区的模糊边缘以及抗噪声影响方面有更好的表现。

图3 细节映射对比

为了更加直观地呈现出细节映射中细节区和非细节区的差异，利用本文的递归阈值选择方法在细节映射上计算阈值TDetail，并定义可视化结果如下：

根据前文分析，非细节区高度相似于真实的云区。基于此，结合待定云区ZInitialn，待定云区精化结果ZFinaln可视化结果如下：

在图4中，给出了结合本文的差异映射、细节映射的整个云检测过程。相比于图4（b）中文献［3］方法得到的待定云区，图4（f）给出的本文方法的结果对于云区的模糊边缘保持得更好。在图4（g）和图4（c）中，本文的细节映射比文献［3］方法得到的细节映射结果，更加精确地抽取了图像的细节。而文献［3］方法得到的细节映射中会含有不应出现的云区信息。图4（d）中，在没有补洞和去除微小聚类前提下，文献［3］方法得到待定云区精化结果显然不能令人满意。即使通过补洞和去除微小聚类等再精化处理后，如图4（e）中所示，相较于真实的云区，文献［3］方法得到云检测结果仍然会出现失真的现象，无法有效保边缘，同时也丢失了一些云区的信息。而在图4（h）中，本文方法的结果很好地反映了真实的云区。

图4 待定云区精化结果对比

2 实验结果与讨论

本文采用2.66 GHz四核CPU、GeForce GTX 480 GPU和4 GB内存的硬件环境，基于CUDA编程实现提出的方法。为了表明提出方法的有效性，将本文方法与文献［3］的彩色航拍图像云检测方法、一些主流的自动图像分割方法［7-9］和交互式图像分割方法［10-11］进行了视觉效果、定量分析等方面的比较，如图5所示。实验测试数据来自于文献［18］中的数据库。

从图5可以看到，文献［3］方法对于复杂的航拍图像不能获得令人满意的结果，而本文方法可以达到很好的效果。

图5 本文方法与文献［3］方法的结果对比

为了客观地评价本文方法的有效性，使用误差率来评估本文方法的精确度。误差率定义为误检测为非云区的云像素的个数以及误检测为云区的非云像素的个数之和与输入的彩色航拍图像像素数量N的比值。在表1中，给出了图5中所涉及结果的误差率。可以看出，本文方法具有更低的误差率。边缘敏感递归滤波是本文方法的核心部件，其时间复杂度是O（N）。借助于GPU的加速计算，本文方法提供毫秒量级的处理反馈。而文献［3］方法含有大量额外的后处理操作，因此，本文方法比文献［3］方法效率更高。同样在表1中，给出图5中所涉及结果的运行时间对比。

表1 图5中结果的误差率和运行时间

图6给出本文方法与自动图像分割方法［7-9］的比较。文献［7］和文献［8］分别具体采用均值漂移和基于图的分割对图像进行分割，然后根据分割区域的平均强度来检测云区。文献［9］方法的结果是基于图像的强度设置两个聚类中心获得。与这些方法相比，本文方法的结果更好。

图6 本文方法与自动图像分割方法的结果对比

在表2中，给出图6中所涉及结果误差率的比较和本文方法的运行时间。可以看出，本文方法具有更低的误差率。这是由于这些自动分割方法可能将大量的有较高强度非云像素误认为云像素，尤其是对半透明云像素检测能力偏弱，因此比本文方法有更高的误差率。

在图7中，进一步将本文方法与文献［10］、文献［11］等交互式图像分割方法比较。由于图7中的航拍图像有大量孤立的小云地区，文献［10］方法、文献［11］方法必须采取大量的人工交互。图7（b）给出了文献［10］方法的人工交互，红色笔画表示云区，蓝色笔画表示非云区。图7（d）给出了文献［11］方法的人工交互，白色笔画表示云区，黑色笔画表示非云区。图7（c）给出了按照图7（b）所示的交互，文献［10］方法的结果。图7（e）给出了按照图7（d）所示的交互，文献［11］方法的结果。本文方法不需要进行交互，仍然获得了满意的结果。

表2 图6中结果的误差率和运行时间

图7 本文方法与交互式图像分割方法的结果对比

表3给出了图7所涉及结果的误差率与交互率的比较和运行时间。其中，交互率为人工交互的笔画所覆盖的像素数量与输入的彩色航拍图像像素数量的比值。可以看出，虽然本文方法不需要交互，仍优于文献［10］方法和文献［11］方法的结果。

表3 图7中结果的交互率、误差率和运行时间

3 结 论

本文提出一种新的边缘敏感的彩色航拍图像云检测框架。首先，计算出差异映射以凸显云区和非云区之间的差异。然后，提出一个新的阈值选择方法，并结合差异映射，获得待定云区。为了有效地从待定云区去除误认为云区的非云像素，将细节映射集成到云检测框架，从而得到最终的云检测结果。为了表明本文方法的有效性，从视觉感受、定量分析等方面比较本文方法和相关工作。实验结果表明，本文方法可以获得令人满意的结果。未来将把学习分类技术应用到云检测中，进一步提高本文方法的准确性。此外，还会将研究工作扩展到多光谱卫星图像的领域上来。

［1］Kang X，Li S，Benediktsson J A.Spectral-spatial hyperspectral image classification with edge-preserving filtering［J］.IEEE Trans.on Geoscience and Remote Sensing，2014，52（5）：2666- 2677.

［2］Frey R A，Ackerman S A，Liu Y，et al.Cloud detection with MODIS.Part I：Improvements in the MODIS cloud mask for collection 5［J］.Journal of Atmospheric and Oceanic Technology，2008，25（7）：1057- 1072.

［3］Zhang Q，Xiao C X.Cloud detection of RGB color aerial photographs by progressive refinement scheme［J］.IEEE Trans.on Geoscience and Remote Sensing，2014，52（11）：7264- 7275.

［4］Farbman Z，Fattal R，Lischinski D，et al.Edge-preserving decompositions for multi-scale tone and detail manipulation［J］.ACM Trans.on Graphics，2008，27（3）：67- 76.

［5］Paris S，Durand F.A fast approximation of the bilateral filter using a signal processing approach［C］∥Proc.of the European Conference on Computer Vision，2006：568- 580.

［6］Yang Q，Tan K H，Ahuja N.Real-time O（1）bilateral filtering［C］∥Proc.of the IEEE Conference on Computer Vision and Pattern Recognition，2009：557- 564.

［7］Comaniciu D，Meer P.Mean shift：a robust approach toward feature space analysis［J］.IEEE Trans.on Pattern Analysis and Machine Intelligence，2002，24（5）：603- 619.

［8］Felzenszwalb P F，Huttenlocher D P.Efficient graph-based image segmentation［J］.International Journal of Computer Vision，2004，59（2）：167- 181.

［9］Lin W T，Lin C H，Wu T H，et al.Image segmentation using the kmeans algorithm for texture features［J］.World Academy of Science，Engineering and Technology，2010，65（4）：612- 615.

［10］Boykov Y Y，Jolly M P.Interactive graph cuts for optimal boundary＆region segmentation of objects in ND images［C］∥Proc.of the International Conference on Computer Vision，2001：105- 112.

［11］Levin A，Lischinski D，Weiss Y.A closed-form solution to natural image matting［J］.IEEE Trans.on Pattern Analysis and Machine Intelligence，2008，30（2）：228- 242.

［12］He K M，Sun J，Tang X O.Guided image filtering［J］.IEEETrans.on Pattern Analysis and Machine Intelligence，2013，35（6）：1397- 1409.

［13］Gastal E S L，Oliveira M M.Domain transform for edge-aware image and video processing［J］.ACM Trans.on Graphics，2011，30（4）：69- 79.

［14］Nehab D，Maximo A，Lima R S，et al.GPU-efficient recursive filtering and summed-area tables［J］.ACM Trans.on Graphics，2011，30（6）：176- 186.

［15］Zhang B，Allebach J P.Adaptive bilateral filter for sharpness enhancement and noise removal［J］.IEEE Trans.On Image Processing，2008，17（5）：664- 678.

［16］Oh J，Kim J，Yuan Z，et al.Adjustment of brightness contrast for fatigue decrease in stereo image production［C］∥Proc.of the IEEE International Symposium on Consumer Electronics，2014：1- 3.

［17］Lai Y K，Rosin P.Efficient circular thresholding［J］.IEEE Trans.on Image Processing，2014，23（3）：992- 1001.

［18］Yao B，Yang X，Zhu SC.Introduction to a large-scale general purpose ground truth database：methodology，annotation tool and benchmarks［C］∥Proc.of the Energy Minimization Methods in Computer Vision and Pattern Recognition，2007：169- 183.

［19］Yang Y R.High dynamic range imaging via tone mapping and multiexposure fusion［D］.Xian：Xidian University，2013.（杨颜如.基于色调映射和多曝光率融合的高动态范围图像成像技术［D］.西安：西安电子科技大学，2013.）

Cloud detection of color aerial images by edge-aware recursive filtering

LIAO Bin，FU Zhong-wang
（School of Computer Science and Information Engineering，Hubei University，Wuhan 430062，China）

In order to detect cloud areas of color aerial images effectively，a cloud detection method of color aerial images by edge-aware recursive filtering is proposed.The variation mapping is constructed using adaptive guided weight based edge-aware recursive filtering.Based on the variation mapping，the input aerial image is divided into pending cloud areas and non-cloud areas by intra-class variance based recursive threshold selection.The original aerial image is progressively filtered using the edge-aware recursive filtering and the hierarchical details are effectively extracted.The detail mapping of the aerial image is constructed and used to refine the pending cloud areas.The performance of the detection results is evaluated by the objective indicator.The experimental results show that the proposed method is effective，easy，and of lower computational complexity compared with other similar methods.

image process；image classification；cloud detection；color aerial images；edge-aware recursive filtering

TP 391.41

A

10.3969／j.issn.1001-506X.2015.12.33

廖 斌（197-9- ），男，副教授，博士，主要研究方向为数字图像处理、计算机图形学。

E-mail：goldgoat＠126.com

付忠旺（199-3- ），男，硕士研究生，主要研究方向为数字图像处理。

E-mail：myexperience＠126.com

1001-506X（2015）12-2879-08

2014- 11- 06；

2015- 05- 24；网络优先出版日期：2015- 07- 06。

网络优先出版地址：http：／／www.cnki.net／kcms／detail／11.2422.TN.20150706.1705.011.html

国家自然科学基金（61300125）资助课题