光学遥感影像云检测方法与流程

本发明涉及光学遥感图像处理技术领域，更具体地说，本发明涉及一种光学遥感影像云检测方法。

背景技术：

以资源三号、高分一号、高分二号为代表，国产光学遥感卫星各项设计指标已逐渐达到国际先进水平，对地观测系统日趋完善，卫星影像的数据量迅速增加，市场化程度逐年提高。然而，并非所有遥感影像均可满足影像信息智能化处理的要求，其中一个很重要的因素就是影像上的云层覆盖。

当前云检测领域中，光谱阈值法是最简单有效的算法，它基于云和地物在可见光波段光谱特性差异，通过亮度阈值实现云与非云目标的区分。常用经验阈值或基于最大类间方差(otsu)原理的自动阈值。此类算法快速有效，但精度略低，会不可避免地对积雪、建筑物、裸地等高亮地物产生误判，且无法定性筛选无云影像。充分利用热红外信息的多光谱综合法可有效改进检测效果，但不适用于光学遥感卫星影像。另一类方法通过分析影像上云和地物纹理特征的差异，提取合适的特征或特征组合，如分形维数、灰度共生矩阵、gabor纹理特征等，区分云和地物。但光学遥感影像上云的种类繁多，不同种类云的特征在各个特征空间内的分布都不集中，使用纹理特征进行精确的云区提取有一定难度。一些改进后算法综合利用影像的辐射和纹理特征，以分类的方式得到云、水、晴空、云影等不同类别，k-means、支持向量机、潜在语义模型等是常用的技术。这些算法在一定程度上提高了检测精度，但需要采用大量经过人工解译且种类不同的含云影像作为样本对分类器进行训练，极其耗时、费力，难以满足海量影像自动化处理的需要。此外，基于同一地区时相相近的两幅或两幅以上影像进行云检测也是一类常见的方法。这类方法将云看作影像中的变化目标，使用变化检测的思想去检测云，常与基于单幅影像的云检测算法联合使用，可有效提高检测精度，但缺点是对数据要求较高，影像本身需具备比较精确的地理信息。

技术实现要素：

针对上述技术中存在的不足之处，本发明提供一种光学遥感影像云检测方法，在较少人工参与的前提下，快速获取精确度较高的影像云掩模和含云量，定性识别无云影像，对全色和多光谱影像均适用。

为了实现根据本发明的这些目的和其它优点，本发明通过以下技术方案实现：

本发明提供一种光学遥感影像云检测方法，其包括步骤：

获取待检测影像的亮度图：将待检测影像中为多光谱的影像转换为单波段亮度影像；

粗估亮度双阈值：根据无云、含云样本影像，计算出对应的最高、最低亮度阈值；

计算精确亮度阈值：分析待检测影像直方图，定性筛选无云影像；对于剩余的含云影像，以粗估的亮度双阈值为限定条件执行基于最大类间方差的计算，得到精确亮度阈值；

云区形态学综合：对阈值分割后的云区执行“腐蚀-条件膨胀-腐蚀”形态学运算，消除似云目标引起的噪声点，填充云缝，优化云区轮廓，输出最终云掩模和含云量。

优选的是，将多光谱影像转换为单波段亮度影像，通过如下公式实现：

p(i,j)＝min[r(i,j),g(i,j),b(i,j)]；

其中，p(i,j)代表转换后的亮度图中位于(i,j)的像元的亮度值，r(i,j),g(i,j),b(i,j)分别代表多光谱图像中位于(i,j)的像元的红、绿、蓝波段的亮度值。

优选的是，计算最高、最低亮度阈值，包括步骤：

人工筛选n1张不含云的卫星影像，逐张统计影像灰度直方图，舍弃直方图位于末端占总数a1％的像素，记录末端截断阈值tend；将所有tend按从高到低顺序排列，舍弃最高的b1％，记录剩余的tend的最大值为高阈值tl-high；

人工勾选n2个云场景，逐张统计场景影像灰度直方图，舍弃直方图位于前端占总数a2％的像素，记录前端截断阈值ssta；将所有ssta按从低到高顺序排列，舍弃最低的b2％，记录剩余的ssta的最小值为低阈值tl-low。

优选的是，n1和n2的取值分别大于100；a1、a2以及a3的取值范围分别为0.01～1；b1、b2以及b3的取值范围分别为1～5。

优选的是，计算精确亮度阈值并定性筛选无云影像，包括步骤：

统计待检测影像的直方图的亮度，将亮度大于tl-high的像素所占比例高于a4％的影像认定为无云影像，其余为含云影像；a4与a1的取值相同；

若是无云影像，则检测结束；若是含云影像，则对直方图中位于高阈值tl-high和低阈值tl-low之间的部分执行基于最大类间方差计算，输出精确阈值tl。

优选的是，对阈值分割后的云区执行“腐蚀-条件膨胀-腐蚀”形态学运算，包括步骤：

对于没有被定性认定为无云的影像，根据亮度阈值tl对影像进行分割，亮度高于阈值tl的部分定义为初始云区；

对所述初始云区，检测面积小于k1的云区，定义为高亮噪声，予以删除，标记为非云；

删除高亮噪声后，对云区执行形尺度为k2的形态学膨胀，在所述形态学膨胀的同时判断新增像元的亮度和在膨胀方向上的亮度梯度，并以此作为膨胀的限定条件；

膨胀处理后，检测面积小于k3的非云区，定义为细小云缝，予以删除，标记为云。

优选的是，所述形态学膨胀满足公式：其中，g为新增像元的亮度，在膨胀方向上的亮度梯度，d是取值范围为0.05～0.25的常数。

本发明至少包括以下有益效果：

本发明提供的光学遥感影像云检测方法，将待检测影像转换成全色后，依次进行粗估亮度双阈值、精确亮度阈值计算以获得精确的亮度阈值对云区进行分割，然后对分割后的云区执行形态学运算，消除似云目标引起的噪声点，填充云缝，优化云区轮廓，输出最终云掩模和含云量；在较少人工参与的前提下，整个检测过程可快速获取精确度较高的影像云掩模和含云量，可定性识别无云影像，检测方法简单有效，对全色和多光谱影像均适用。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为本发明所述的光学遥感影像云检测方法的示意图；

图2为本发明所述的光学遥感影像云检测方法的流程图；

图3(a)-图3(b)为粗估最高亮度阈值的示意图；

图4(a)-图4(b)为粗估最低亮度阈值的示意图；

图5(a)-图5(b)为以阈值为限定条件计算精确亮度阈值的示意图；

图6(a)-图6(e)为阈值分割及云区形态学运算过程示意图；

图7(a)-图7(e)为图6(a)-图6(e)一一对应的局部放大示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不排除一个或多个其它元件或其组合的存在或添加。

如图1和图2所示，本发明提供一种光学遥感影像云检测方法，其包括步骤：

s10，获取待检测影像的亮度图：将待检测影像中为多光谱的影像转换为单波段亮度影像；

s20，粗估亮度双阈值：根据无云、含云样本影像，计算出对应的最高、最低亮度阈值；

s30，计算精确亮度阈值：分析待检测影像直方图，定性筛选无云影像；对于剩余的含云影像，以粗估的亮度双阈值为限定条件执行基于最大类间方差的计算，得到精确亮度阈值；

s40，云区形态学综合：对阈值分割后的云区执行“腐蚀-条件膨胀-腐蚀”形态学运算，消除似云目标引起的噪声点，填充云缝，优化云区轮廓，输出最终云掩模和含云量。

上述实施方式中，光学遥感卫星影像一般包含全色影像和多光谱影像(通常有蓝、绿、红和红外4个波段)。步骤s10中，获取待检测影像的亮度图时，先判断待检测影像是全色还是多光谱影像：若是全色，则直接跳至步骤s20，若是多光谱影像，需将待检测影像中为多光谱的影像转换为单波段亮度影像；转换通过如下公式实现：p(i,j)＝min[r(i,j),g(i,j),b(i,j)](1)；其中，p(i,j)代表转换后的亮度图中位于(i,j)的像元的亮度值，r(i,j),g(i,j),b(i,j)分别代表多光谱图像中位于(i,j)的像元的红、绿、蓝波段的亮度值。由于云对太阳光呈现米式散射，对各波段均有比较强的散射，表现在多光谱影像上，云是白色的，各波段亮度值都很高。而非云地表目标对太阳光呈漫反射，不同波段反射率往往不同，呈现在多光谱图像上，非云目标往往是彩色的，各波段亮度值高低不同，但最小亮度值往往较低。因此，通过上述公式(1)提取各波段亮度值的最小值得到的单波段亮度图像，综合了多光谱图像的亮度和饱和度信息，更容易实现云与非云目标的区分。

步骤s20中，通过一定数目的无云样本影像计算最高亮度阈值，是为了保证云的准确率；通过一定数目的云样本影像计算最低亮度阈值，是为了保证云的查全率。作为一种实施例，计算最高、最低亮度阈值，包括步骤：

s21，人工筛选n1张不含云的卫星影像，逐张统计影像灰度直方图，舍弃直方图位于末端占总数a1％的像素，记录末端截断阈值tend；将所有tend按从高到低顺序排列，舍弃最高的b1％，记录剩余的tend的最大值为高阈值tl-high；

s22，人工勾选n2个云场景，逐张统计场景影像灰度直方图，舍弃直方图位于前端占总数a2％的像素，记录前端截断阈值ssta；将所有ssta按从低到高顺序排列，舍弃最低的b2％，记录剩余的ssta的最小值为低阈值tl-low。

其中，n1张不含云的卫星影像，应尽量选取包含多种场景的影像，如植被、城镇、裸地和水面等，但不能选取含有积雪的影像。n2个云场景，应尽量包含薄云和厚云，但必须避开云缝和透光薄雾。n1和n2的取值分别大于100；a1、a2以及a3的取值范围分别为0.01～1；b1、b2以及b3的取值范围分别为1～5。作为一种优选实施例，a1＝a2＝0.1；b1＝b2＝2。对于某一类特定传感器，在前期传感器校正与相对辐射定标工作没有较大改动情况下，步骤s21和步骤s22中涉及的两阈值tl-high和tl-low具有普适性，即适用于所有同类型传感器影像。图3(a)至图3(b)给出了以高分一号全色影像为例的粗估最高亮度阈值示意图，其中，图3(a)为影像缩略图示意，图3(b)为影像直方图及末端截断阈值示意。图4(a)至图4(b)给出了以高分一号全色影像为例的粗估最低亮度阈值示意图，其中，图4(a)为影像缩略图示意；图4(b)为影像直方图及前端截断阈值示意。

在步骤s20粗估亮度双阈值的基础上，步骤s30用于进一步精确计算亮度阈值。作为一种实施例，计算精确亮度阈值并定性筛选无云影像，包括步骤：

s31，统计待检测影像的直方图的亮度，将亮度大于tl-high的像素所占比例高于a4％的影像认定为无云影像，其余为含云影像；a4与a1的取值相同。

s32，若是无云影像，则检测结束；若是含云影像，则对直方图中位于高阈值tl-high和低阈值tl-low之间的部分执行基于最大类间方差计算，输出精确亮度阈值tl。

图5(a)至图5(b)给出了以高分一号全色影像为例的、双阈值为限定条件计算精确亮度阈值示意图；其中，图5(a)为影像缩略图示意，图5(b)为影像直方图及灰度阈值示意，两侧直线代表灰度限定条件，中间直线代表带限定条件的精确亮度阈值。

作为另一种实施例，对阈值分割后的云区执行“腐蚀-条件膨胀-腐蚀”形态学运算，包括步骤：

s33，对于阈值分割后的初始结果，检测面积小于k1的云区，定义为高亮噪声，予以删除，标记为非云；

s34，删除高亮噪声后，对云区执行形尺度为k2的形态学膨胀，在形态学膨胀的同时判断新增像元的亮度和在膨胀方向上的亮度梯度，并以此作为膨胀的限定条件；

s35，膨胀处理后，检测面积小于k3的非云区，定义为细小云缝，予以删除，标记为云。

上述步骤s34中的形态学膨胀满足公式：其中，g为新增像元的亮度，在膨胀方向上的亮度梯度，d是取值范围为0.05～0.25的常数。作为进一步优选，d＝0.15。若满足公式(2)列出的限定条件，才执行膨胀，否则不膨胀。步骤s33、s34以及s35中涉及的k1、k2和k3均为配置参数，实际应用中可根据具体情况改动，例如若场景中常含有大型军事目标(机场、靶场等)，应设定一个较大的k1来避免错检；若需要充分挖掘云缝中的有效信息，应设置一个较小的k2和k3；若重点关注云区的查全率，不希望得到过于破碎的云掩模，应设置一个较大的k2和k3。图6(a)至图6(e)给出了阈值分割及云区形态学运算过程的示意，图7(a)至图7(e)为图6(a)至图6(e)一一对应的局部放大示意图。其中，图6(a)为影像缩略图示意，图6(b)为阈值分割后的云区示意，图6(c)为去除小面积噪声点后的云区示意，图6(d)为带限定条件的形态学膨胀后的云区示意，图6(e)为填充掉小面积云层缝隙后的云区示意。

本发明提供的光学遥感影像云检测方法，将待检测影像转换成全色后，依次进行粗估亮度双阈值、精确亮度阈值计算以获得精确的亮度阈值对云区进行分割，然后对分割后的云区执行形态学运算，消除似云目标引起的噪声点，填充云缝，优化云区轮廓，输出最终云掩模和含云量；整个检测过程可快速获取精确度较高的影像云掩模和含云量，可定性识别无云影像，检测方法简单有效，对全色和多光谱影像均适用。对于某一指定传感器影像，在较少人工参与的前提下，即可实现海量影像快速、自动化检测，可满足实际生产的需要。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。