一种快速的遥感卫星图像云检测算法∗

关峥嵘 谭毅华 田金文

（华中科技大学自动化学院 武汉 430074）

1 引言

随着遥感技术的不断发展，我们能获得的遥感卫星影像数据量也越来越多，遥感卫星影像的分辨率也越来越高，这些卫星影像携带的地面目标信息也越来越多。但同时，地球表面是包裹着大气层，遥感传感器常常受到大气或者云层的干扰［1］，使得我们对地球表面的观测受到遮挡，因此针对遥感卫星影像的云检测也就成为了遥感卫星图像处理的关键性技术之一，云检测也被认为是遥感影像后续处理（识别和分类等）的预处理技术。此外在气象学方面，针对云的检测也被认为是进行气象预测的重要手段之一［2］。

目前来说，针对云检测是方法大致分为三种：阈值法，纹理分析法和基于机器学习的方法。基于阈值的方法利用了不同物质反应在图像上不同的明暗程度，又可以分为简单的单阈值［3～4］，动态阈值法［5］；纹理分析法从本质上来说就是利用同类物质之间的相似性和不同物质之间的不连续性，往往需要结合其他信息一起进行判别，比如多时相［6～7］和物体反射率［8］，这些方法虽然都取得了一定的进步，但是大都适应性不强，利用了一些先验知识，导致了在场景复杂的情况下的检测效果较差。而基于机器学习的方法往往能获得比较好的效果，利用训练样本来进行特征训练，使得算法的适应性大大加强，支持向量机（Support Vector Machine，SVM）是最常用的分类器［9～12］，还有利用神经网络模型来进行训练［13］。此外一些比较新颖的图像处理技术也被应用到云检测当中，比如超像素分割［14～15］，模糊聚类［16～17］。利用机器学习来进行云检测时，检测效果往往取决于选取的特征描述能力的强弱，因此如何寻找出一种好的描述特征，是我们目前面临的挑战之一。论文针对该问题，提出了一种基于改进MR 的描述特征来进行云检测，并在高分一号卫星图像数据集上获得了令人满意的效果。

2 最大响应滤波器（Maximum Response Filter，MR）的快速实现

2.1 MR的基本原理

最大响应滤波器（Maximum Response Filter）［18］采用了38个滤波器，但最终只选取8组最大的响应值作为输出，因此又被称为MR8。MR8 的38 个滤波器包含了18 个高斯一阶导数滤波器（三个尺度六个方向），18 个高斯二阶倒数滤波器（三个尺度六个方向），尺度为(σx，σy)={(1 ，3)，( 2，6)，(4，12)}，方向角取值为θ={0°，30°，60°，90°，120°，150°}；1个高斯滤波器和1个高斯型的拉普拉斯滤波器（Laplacian of Gaussian），标准差σ 取10。38 个滤波器的时域图如图1所示。

由于MR8 集合了高斯一阶和二阶导数滤波器，因此能有效地提取图像的边缘信息，此外由于采用了不同尺度的滤波器使得提取的特征能适应不同的图像分辨率，不同的方向滤波器使得其在非均质的纹理下获得更强的分类能力，而最终只选取最大的8个响应值，使得MR8特征具有旋转不变的性质。综上所述，MR8特征具有较强的特征描述能力。

2.2 MR8快速实现方法

在计算MR8 特征的过程中，由于需要使用高斯核与图像进行大量的卷积操作，而通常遥感卫星图像分辨率较高，幅面较大，导致整个提取特征的过程中非常耗时，因此如何快速地实现MR8 是急需解决的问题。受之前的工作启发［19～20］，本论文使用递归（Infinite Impulse Response，IIR）数字滤波器来快速近似，以设计一种快速实现MR8的方法。

一维信号的高斯核卷积过程可用两次递归的形式来实现，第一个是前向递归（Forward，n递增）：

第二个是后向递归（Backward，n递减）：

可以表示为如图2所示。

图2 一维信号的卷积通过递归来实现

式（1）和式（2）中的参数B，b0，b1，b2，b3，由下面的式子决定：

其中m0，m1，m2为常数，取值分别为1.16680，1.10783，1.40586。而q 是与标准差σ 有关的数值，一般按照如下方法取值：

根据上面的方法把一维信号的高斯卷积操作变成了两次递归来实现，这很容易推广到二维图像的高斯卷积操作上去，最后只需要做四次递归操作即可：从上到下，从左到右，从下到上和从右到左。

2.3 图像块特征选取

为了尽可能好地描述出云的图像特征，我们这里并没有直接采用原始的MR8 特征进行操作，而是利用了K-D 树来构建MR8 直方图特征用于判别。具体的步骤如下所示。

1）计算训练数据集上的MR8 特征，并利用K均值聚类算法得到纹理特征的数据字典；

2）利用得到的数据字典构建K-D树；

3）得到图像块的MR8 特征，将得到的纹理特征进行K-D树查询；

4）生成查询结果的归一化直方图，并利用该直方图特征作为云判别的特征。

2.4 算法流程设计

本文提出的算法步骤可以如图3所示：

图3 算法流程图

在训练阶段，分为两步：第一步是构建数据字典，首先是提取图像块（大小为128*128）的MR8特征，进行K均值聚类以后为了避免特征维数过大导致的维数灾难，我们这里使用了主成分分析（Principal Component Analysis，PCA）进行了降维处理，最后得到的这些聚类中心就构成了纹理特征的数据字典，然后利用数据字典构建K-D 树，如果按照传统的方法在得到数据字典后需要与聚类中心进行欧式距离遍历，再去统计得到结果的直方图特征，这样的操作不仅耗时而且特征的描述效果也不甚理想；第二步是利用数据字典得到图像块的纹理特征，随后把这些纹理特征利用支持向量机（Support Vector Machine，SVM）来训练出一个分类器。

在测试阶段时与上面类似，只需要得到图像块的直方图特征，然后利用SVM去判别即可。

3 实验结果

3.1 实验数据准备

在实验过程中，训练所采用的数据来自与谷歌地图（Google Map），测试所用的数据来自于高分一号卫星数据（拍摄于云南某地，分辨率为2m）。训练数据集一共有6000 个样本，包括2000 个正样本（云区图像），4000 个负样本（包括森林，湖泊，城市，草地等区域图像）。

为了确定实验中的数据块大小、词典的单词数量和是否采取PCA等参数，我们做了一系列的实验来进行对比。

根据实验结果，最终我们选定使用128*128 的图像块大小，数据字典数量为1024，采用51*51 的高斯卷积核，因此我们在图像上得到的MR8 特征其实是78*78*8 维的特征，再进行K-D 树查询后，得到了78*78 维的查询结果。接下来对其进行量化并统计其直方图分布，作为该图像块的描述特征。从实验中我们可以看到，进行PCA操作有效地缩短了算法所需的时间，而且去除了冗余数据后也使得检测精度有小范围的提升。

图4 高分一号卫星的测试遥感图像示例

表1 不同的单词个数对准确率的影响（图像块大小为128*128）

表2 不同的图像块对准确率的影响（单词数量为1024）

表3 PCA对时间和准确率的影响（2048*2048大小图像）

3.2 实验结果对比

测试图像我们选用了35 幅2048*2048 大小的卫星图像，并在不同的算法上进行了对比，其中我们用F_alarm 表示虚警率，准确率precision 通过TP/(TP+FP)得到，召回率通过TP/(TP+FN)得到，精度accuracy 通过(TP+TN)/(TP+FN+TN+FP)计算得到。

表4 不同方法的云检测结果对比

从表中可以看到，我们的算法以微小的时间代价在各方面都得到了最好的结果，CI 和DT 的方法虽然耗费时间比较短，但是损失了较大的准确率，因此可以认为本论文提出的算法是优于这些方法的。

4 结语

本文提出了一种通过快速MR8 和K-D 树构建直方图特征来进行云检测的算法，与一般的检测方法不同的是，该方法通过四个方向上的递归操作来快速实现MR8 特征的计算，并利用数据字典和K-D 树构建直方图特征，基于MR8 的直方图特征能有效地描述云块区域，具有较强的适应性而且执行速度较快，能满足不同条件下快速进行遥感卫星图像的自动云判别的要求。