本发明属于遥感图像处理技术领域，尤其涉及一种基于空域自然场景统计的遥感图像云检测方法。
背景技术：
遥感是上世纪60年代发展起来的与电子、空间、计算机、光学、地理等科学技术密切相关的一门新兴的交叉科学技术。它不仅是当前信息领域中一个十分活跃的学科，也是当今科学技术的一个重要组成部分。以资源三号、高分一号、高分二号为代表，国产光学遥感卫星各项设计指标已逐渐达到国际先进水平，对地观测系统日趋完善，卫星影像的数据量迅速增加，市场化程度逐年提高。然而，并非所有遥感影像均可满足影像信息智能化处理的要求，其中一个很重要的因素就是影像上的云层覆盖。当有大量云层覆盖的情况下，卫星成像设备无法正确接收到来自地面目标的光反射特性，相应的可见光卫星图像数据包含的地面目标可用信息大幅减少，对后续基于图像分类、目标检测及识别等处理和应用产生了极大的负面影响，降低了可见光卫星遥感图像的使用价值。综上所述，云检测是可见光卫星图像处理应用的首要任务及关键技术之一。
现在的云检测算法可以分为两类：不同时间的云检测算法和同一时间的云检测算法。不同时间的云检测，这类算法主要遵循的原则是在短时间内遥感地物的反射率和温度变化不大，2013年，sumingjin提出使用两天内的同一地点的图像，来对比他们之前的差异，从而检测云的位置，nicholasr.goodwin通过提取不同气候图像的光谱和纹理特征信息来和待检测图像进行对比，从而实现图像的云检测算法。2014年，zhezhu提出要使用多幅不同时间的遥感图像，建立了一个多时间滤波器，从而实现遥感图像云区域的自动检测。然而上面提到的方法，计算复杂度高，检测结果严重依赖于配准的效果，很难实际应用于卫星上面。同一时间的云检测算法，这类算法主要遵循的原理是云通常是比地面更亮和温度更低的，这类算法又可以分为两种，一种是基于多通道图像特性的算法，另一种是基于可见光通道图像特性的算法。多通道图像特性的算法，需要使用中波红外通道来设置温度阈值，来尽可能的区分云与其他区域，但是这就限制了该类算法的应用环境，例如中国发射的卫星资源三号、高分一号、高分二号都没有中波红外波段，所以上述算法都失效了。另外一种算法需要，提取可见光图像信息，例如h.k.chethan于2009年提出使用gabor滤波器来提取云的纹理信息，并通过svm对云的图像进行分类。在2015年，yiyuan提出使用超像素来对图像进行分割，并将分割的结果通过svm来对图像进行分类，得到最终的云检测结果。云检测算法方法众多，但是目前仍然有以下几个难点需要着重解决，
1)碎云形状多样，对基于纹理和结构的云检测方法带来了挑战；
2)云有薄厚的区分，简单的阈值分割会造成薄云的检测丢失；
3)相似地物的干扰，如高亮建筑、雪等都会给云检测带来影响。
技术实现要素：
为解决上述问题，本发明提供一种基于空域自然场景统计的遥感图像云检测方法，能够准确地将云区域从可见光波段遥感图像中提取出来。
一种基于空域自然场景统计的遥感图像云检测方法，包括以下步骤：
s1：将可见光波段遥感图像进行超像素分割，得到超像素块；
s2：分别获取各超像素块对应的非对称广义高斯分布的波形参数f1与云图像对应的非对称广义高斯分布的波形参数f2、无云图像对应的非对称广义高斯分布的波形参数f3之间的距离，如果波形参数f1与波形参数f2之间的距离，小于波形参数f1与波形参数f3之间的距离，则波形参数f1对应的超像素块为云雪区域，否则为非云雪区域；
s3：获取所述云雪区域对应的超像素块的gabor特征；
s4：将所述gabor特征输入支持向量机svm中进行判断，完成对云区域的检测；其中，所述支持向量机svm根据云图像的gabor特征与雪图像的gabor特征训练得到。
可选的，所述将可见光波段遥感图像进行超像素分割具体为：
采用简单线性迭代聚类法对可见光波段遥感图像进行超像素分割，其中，在简单线性迭代聚类法的增强连通性的过程中，如果面积小于设定阈值t2的超像素块的像素均值大于设定阈值t3，则面积小于设定阈值t2的超像素块不并入与其最近的面积不小于设定阈值t2的超像素块；如果面积小于设定阈值t2的超像素块的像素均值不大于设定阈值t3，则面积小于设定阈值t2的超像素块并入与其最近的面积不小于设定阈值t2的超像素块。
可选的，在将可见光波段遥感图像进行超像素分割，获取超像素块后，剔除面积小于设定阈值t4的超像素块，再将剩余的超像素块执行步骤s2和步骤s3。
可选的，所述获取各超像素块对应的非对称广义高斯分布的波形参数f1具体为：
获取超像素块各像素点的像素值、局部均值以及局部标准差；
根据各像素点的像素值、局部均值以及局部标准差，获取各像素点的均差对比度归一化系数；
根据各像素点的均差对比度归一化系数，获取各超像素块对应的均差对比度归一化系数矩阵；
根据所述均差对比度归一化系数矩阵，获取各超像素块对应的非对称广义高斯分布，并获取非对称广义高斯分布的波形参数f1。
可选的，所述波形参数f1与波形参数f2之间的距离的获取方法为：
获取超像素块对应的非对称广义高斯分布的波形参数，其中，所述波形参数包括波形半高宽、左偏移量、右偏移量以及均值；
将超像素块各像素点作为中心像素点，则将各中心像素点的均差对比度归一化系数分别与中心像素点的水平方向、竖直方向、对角方向以及反角方向的像素点的均差对比度归一化系数相乘，从而获取超像素块在水平方向、竖直方向、对角方向以及反角方向对应的均差对比度归一化系数矩阵；
根据超像素块在水平方向、竖直方向、对角方向以及反角方向对应的均差对比度归一化系数矩阵，分别获取超像素块在水平方向、竖直方向、对角方向以及反角方向对应的非对称广义高斯分布；
分别获取超像素块在水平方向、竖直方向、对角方向以及反角方向对应的非对称广义高斯分布的波形参数；
获取超像素块20个特征值对应的协方差矩阵；其中，所述20个特征值分别为超像素块对应的非对称广义高斯分布的波形参数，以及超像素块在水平方向、竖直方向、对角方向、反角方向对应的非对称广义高斯分布的波形参数；
根据超像素块的20个特征值、超像素块的20个特征值对应的协方差矩阵、云图像的20个特征值、云图像的20个特征值对应的协方差矩阵，获取超像素块对应的非对称广义高斯分布的波形参数f1与云图像对应的非对称广义高斯分布的波形参数f2之间的距离。
可选的，所述波形参数f1与波形参数f3之间的距离的获取方法为：
根据超像素块的20个特征值、20个特征值对应的协方差矩阵、非云图像的20个特征值、20个特征值对应的协方差矩阵，获取超像素块对应的非对称广义高斯分布的波形参数f1与非云图像对应的非对称广义高斯分布的波形参数f3之间的距离。
可选的，所述获取所述云雪区域对应的超像素块的gabor特征具体为：
在0°～３６0°的范围中随机选取两个以上的角度作为gabor滤波器的核函数的方向；
根据经验值为所述核函数确定两个以上的正弦函数波长，其中所述正弦函数波长确定所述核函数的频率；
根据所述方向和频率，获取所述云雪区域对应的超像素块在不同方向和频率上的gabor值；
获取多个gabor值的平均值，将所述平均值作为所述云雪区域对应的超像素块的gabor特征。
有益效果：
1、本发明提供一种基于空域自然场景统计的遥感图像云检测方法，由于云区域比较平坦，非对称广义高斯分布的值集中在零附近，所以云区域的非对称广义高斯分布近似尖峰状态；而城市高亮度区域，比较分散，所以非对称广义高斯分布比较平缓；则当待测超像素块的非对称广义高斯分布与云图像的非对称广义高斯分布一致，即波形参数f1与波形参数f2之间的距离，小于波形参数f1与波形参数f3之间的距离，那么待测超像素块判断为云雪区域，否则判断为非云雪区域；这种统计参数距离的方法首次应用于云检测，能够准确地将云区域从可见光波段遥感图像中提取出来；
由于云是漂浮在空中，由汽化的水汽不断聚合而形成的，则云周围纹理过渡平滑，而雪是紧贴地表周围过渡明显，纹理较大；同时汽化的云的纹理方向没有规律，而由于受地面纹理的影响，雪区域纹理有明显的方向性。基于这个两个性质，本发明提出使用gabor特征来进行云雪区域的区分，gabor特征通过gabor滤波器获得，gabor滤波器又能很好的提取特定方向的纹理信息，因此可以用来区分云雪边缘的纹理强弱；最后本发明再通过svm来区别云区域和雪区域，使得本发明的对云雪的区别明显，云检测率高，满足卫星系统在轨云评估的问题。。
2、本发明使用简单线性迭代聚类(slic)来对可见光波段遥感图像进行超像素分割，传统的slic为了保证分割质量，会把小面积区域合并到临近区域，这样就会把碎云合并到其它区域；而本发明在slic增强连通性的过程中，对小面积区域是否合并入大面积的临近区域进行判断，如果小面积的超像素块的像素均值大于设好的阈值，则高亮度的小面积超像素块不合并到其他区域，这样就保留了碎云，提高了分割的结果准确性。
3、本发明提取不同方向和频率的gabor值，将多个gabor值的平均值作为gabor特征，这样使得云雪区域在不同方向上有所区分，能够进一步提高云检测率。
附图说明
图1为本发明提供的一种基于空域自然场景统计的遥感图像云检测方法的流程图；
图2为本发明提供的一种基于空域自然场景统计的遥感图像云检测方法的流程效果图；
图3为本发明提供的云图像、无云图像的均差对比度归一化系数分布曲线；
图4为本发明提供的中心像素点与其水平、竖直、对角以及反角方向的像素点的示意图；
图5(a)为本发明提供的可见光波段云区域遥感图像示意图；
图5(b)为本发明提供的可见光波段雪区域遥感图像示意图；
图6(a)为本发明提供的可见光波段云区域遥感图像示意图对应的梯度统计直方图；
图6(b)为本发明提供的可见光波段雪区域遥感图像示意图对应的梯度统计直方图；
图7(a)为本发明提供的可见光波段遥感图像示意图；
图7(b)为本发明提供的可见光波段遥感图像示意图对应的gabor特征示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。
参见图1，该图为本实施例提供的一种基于空域自然场景统计的遥感图像云检测方法的流程图。图2为本实施例提供的一种基于空域自然场景统计的遥感图像云检测方法的流程效果图。
一种基于空域自然场景统计的遥感图像云检测方法，包括以下步骤：
s1：将可见光波段遥感图像进行超像素分割，得到超像素块。
下面介绍采用简单线性迭代聚类法slic将可见光波段遥感图像进行超像素分割的方法，包括以下步骤：
s101：初始化种子点(聚类中心)：按照设定的超像素块个数，在可见光波段遥感图像内均匀的分配种子点。假设可见光波段遥感图像总共有n个像素点，预分割为k个相同尺寸的超像素块，那么每个超像素块的大小为n/k，则相邻种子点的距离(步长)近似为s＝sqrt(n/k)(自动整除)。
s102：在种子点的n×n邻域内重新选择种子点(一般取n＝3)。具体方法为：计算该邻域内所有像素点的梯度值，将种子点移到该邻域内梯度值最小的地方。这样做的目的是为了避免种子点落在梯度较大的轮廓边界上，以免影响后续聚类效果。
s103：在每个种子点周围的邻域内为每个像素点初始化分配类标签(即属于哪个聚类中心)。与标准的k－means聚类方法在整张可见光波段遥感图像中搜索不同，本实施例采用的简单线性迭代聚类法slic的搜索范围限制为2s×2s，可以加速算法收敛。
s104：距离度量。包括颜色距离和空间距离，对于每个搜索到的像素点，分别计算它和该种子点的距离。距离计算方法如下：
其中，可见光波段遥感图像u，v的颜色信息为(lu，au，bu)，(lv，av，bv)，lab颜色空间是由rgb空间线性变换后的空间，这样的处理有利于颜色信息的提取，空间位置信息为(xu，yu)，(xv，yv)，dc代表颜色距离，ds代表空间距离，ns是类内最大空间距离，定义为ns＝s＝sqrt(n/k)，适用于每个聚类。最大的颜色距离nc，随可见光波段遥感图像不同而不同，也随聚类不同而不同，所以本发明取一个固定常数m(取值范围[1，40]，一般取10)代替。最终的距离度量d＇如下：
由于每个像素点都会被多个种子点搜索到，所以每个像素点都会有一个与周围种子点的距离，取最小值对应的种子点作为该像素点的聚类中心。
s105：迭代优化。理论上上述步骤不断迭代直到误差收敛(可以理解为每个像素点聚类中心不再发生变化为止)，实践发现10次迭代对绝大部分可见光波段遥感图像都可以得到较理想效果，所以一般迭代次数取10。
s106：增强连通性。经过上述迭代优化可能出现以下瑕疵：出现多连通情况、超像素块尺寸过小，单个超像素块被切割成多个不连续超像素等，这些情况可以通过增强连通性解决，判断这个区域大小是否与周围区域连同，如果不连通，直接合并到其他区域上，在这个合并的基础上，本发明添加如下判断，如果面积小于设定阈值t2的超像素块的像素均值大于设定阈值t3，则面积小于设定阈值t2的超像素块不并入与其最近的面积不小于设定阈值t2的超像素块；如果面积小于设定阈值t2的超像素块的像素均值不大于设定阈值t3，则面积小于设定阈值t2的超像素块并入与其最近的面积不小于设定阈值t2的超像素块，从而保留对碎云的检测。
可选地，在将可见光波段遥感图像进行超像素分割，获取大小不一的超像素块后，剔除小于设定阈值t4的超像素块，再将剩余的超像素块执行步骤s2和步骤s3。由此可见，剩余的超像素块面积大小不尽相同。
可选地，采用大津阈值分割获取设定阈值t4，包括以下步骤：
其中，t为大津阈值，
s2：分别获取各超像素块对应的非对称广义高斯分布的波形参数f1与云图像对应的非对称广义高斯分布的波形参数f2、无云图像对应的非对称广义高斯分布的波形参数f3之间的距离，如果波形参数f1与波形参数f2之间的距离，小于波形参数f1与波形参数f3之间的距离，则波形参数f1对应的超像素块为云雪区域，否则为非云雪区域。
需要说明的是，云图像和无云图像可以由已知图像库提供，且该已知图像库提供的云图像和无云图像为经过标定的图像。
可选地，所述各超像素块对应的非对称广义高斯分布的获取方法为：
s201：获取超像素块各像素点的像素值、局部均值以及局部标准差。
具体地，各像素点的局部均值μ(i,j)以及局部标准差σ(i,j)计算方法分别为：
其中，i,j是像素点的位置，如果超像素块的大小为m0×n0，则有i∈{1,2,…,m0}，j∈{1,2,…,n0}，ω＝{ωk,l|k＝-k,…,k,l＝-l,…,l}是一个二维对称的单位高斯权重核，一般设置k＝l，取值根据超像素块大小而定。常用的设置为k＝l＝3。
s202：根据各像素点的像素值、局部均值以及局部标准差，获取各像素点的均差对比度归一化系数。
具体地，各像素点的均差对比度归一化系数mscn的计算方法为：
其中，i(i,j)表示归一化前像素点的像素值，
参见图3，该图为本实施例提供的云图像、无云图像的均差对比度归一化系数分布曲线。由于云雪区域比较平坦，均差对比度归一化后，区域的值集中在零附近，所以云雪区域的均差对比度归一化系数分布近似尖峰状态，如图3中虚线，城市高亮度区域，即非云雪区域比较分散，均差对比度归一化系数分布近似高斯状态，如图3中实线圆圈分布。由此可知，云雪区域对应的分布曲线较高，非云雪区域对应的分布曲线低矮，数值比较分散。当待测超像素块的均差对比度归一化系数分布与云雪区域的分布一致，我们就判断为云雪区域的超像素块，否则判断为非云雪超像素块。
s203：根据各像素点的均差对比度归一化系数，获取各超像素块对应的均差对比度归一化系数矩阵。
s204：根据所述均差对比度归一化系数矩阵，获取各超像素块对应的非对称广义高斯分布。
具体地，非对称广义高斯分布
其中，γ是非对称广义高斯分布的波形半高宽(波形的“胖瘦”)，βl,βr分别是波形左偏移量、右偏移量，可以通过矩匹配的方法求解得到。同时，分布的均值η也对非对称广义高斯分布有影响，因此也作为非对称广义高斯分布建模的特征之一：
由图3所示的分布曲线可知云雪区域的均差对比度归一化系数分布不是一个规则的高斯分布。为了计算分布的偏差，需要对水平、竖直、对角以及反角方向的临近系数乘积，再将乘积的结果作非对称广义高斯分布的建模，这样就可以将超像素块像素间的结构关系进行表征。参见图4，该图为本实施例提供的中心像素点与其水平、竖直、对角以及反角方向的像素点的示意图。由图4可知，中心像素点与其水平、竖直、对角以及反角方向的临近系数乘积分别为：
可选地，所述波形参数f1与波形参数f2，波形参数f1与波形参数f3之间的距离的获取方法为：
s205：获取超像素块对应的非对称广义高斯分布的波形参数，其中，所述波形参数包括波形半高宽、左偏移量、右偏移量以及均值。
s206：将超像素块各像素点作为中心像素点，则将各中心像素点的均差对比度归一化系数分别与中心像素点的水平方向、竖直方向、对角方向以及反角方向的像素点的均差对比度归一化系数相乘，从而获取超像素块在水平方向、竖直方向、对角方向以及反角方向对应的均差对比度归一化系数矩阵。
s207：根据超像素块在水平方向、竖直方向、对角方向以及反角方向对应的均差对比度归一化系数矩阵，分别获取超像素块在水平方向、竖直方向、对角方向以及反角方向对应的非对称广义高斯分布。
s208：分别获取超像素块在水平方向、竖直方向、对角方向以及反角方向对应的非对称广义高斯分布的波形半高宽、左偏移量、右偏移量以及均值。
s209：获取超像素块20个特征值对应的协方差矩阵；其中，所述20个特征值分别为超像素块、以及超像素块在水平方向、竖直方向、对角方向、反角方向对应的非对称广义高斯分布的波形半高宽、左偏移量、右偏移量以及均值。
s210：根据超像素块的20个特征值、超像素块的20个特征值对应的协方差矩阵、云图像的20个特征值、云图像的20个特征值对应的协方差矩阵，获取超像素块对应的非对称广义高斯分布的波形参数f1与云图像对应的非对称广义高斯分布的波形参数f2之间的距离。
s211：根据超像素块的20个特征值、20个特征值对应的协方差矩阵、非云图像的20个特征值、20个特征值对应的协方差矩阵，获取超像素块对应的非对称广义高斯分布的波形参数f1与非云图像对应的非对称广义高斯分布的波形参数f3之间的距离。
例如，基于上述距离获取方法，先在已知图像库中求出云雪区域，即云图像的特征v1(20个特征)、云图像特征的协方差矩阵∑1，非云雪区域，即非云图像的特征v2、非云图像特征的协方差矩阵∑2，再求待检测的超像素块的特征v3，超像素块特征的协方差矩阵∑3，带入到下面的求距离公式，如果d1<d2，其中，d1为波形参数f1与波形参数f2之间的距离，d2为波形参数f1与波形参数f3之间的距离，则判断待检测的超像素块为云雪区域，否则为非云雪区域。
s3：获取所述云雪区域对应的超像素块的gabor特征。
需要说明的是，gabor特征是由gabor滤波器产生的，可见光波段遥感图像经过gabor滤波器后会产生一个gabor值，该值即为gabor特征。gabor滤波器是一类十分重要的方向性滤波器，特别适合在二维图像纹理分析中检测方向性较强的边缘。gabor滤波器具有方向和频域可调的性质，其核函数g(x,y；λ,θ,ψ,σ,ε)如下所示：
x′＝xcosθ+ysinθ
y′＝-xsinθ+ycosθ
其中λ是正弦函数波长，决定gabor核函数的频率，θ是gabor核函数的方向，ψ是相位偏移，σ高斯函数的标准差，ε是空间的宽高比。由此可见，通过将云雪区域对应的超像素块通过gabor滤波器，并为gabor滤波器核函数设置不同的方向与正弦函数波长，则能够获取云雪区域对应的超像素块在不同方向和频率上的gabor值。然后获取多个gabor值的平均值，将所述平均值作为所述云雪区域对应的超像素块的gabor特征。
可选地，在本实施例中，gabor核函数的方向θ与正弦函数波长λ具体为：
θ∈{0°,45°,90°,135°,180°,225°,270°,315°},λ∈{2,4}
需要说明的是，gabor核函数的方向除了采用上述设置方法，还可以在0°～360°的范围中随机选取两个以上的角度，如4个角度、6个角度或者9个角度等。同理，gabor核函数的正弦函数波长也可以根据实际需要或者根据经验值选取，本实施例对此不作赘述。
s4：将所述gabor特征输入支持向量机svm中进行判断，完成对云区域的检测；其中，所述支持向量机svm根据云图像的gabor特征与雪图像的gabor特征训练得到。
需要说明的是，支持向量机svm是建立在计学习理论的vc维理论和结构风险最小原理基础上的，根据有限的样本信息在模型的复杂性(即对特定训练样本的学习精度)和学习能力(即无错误地识别任意样本的能力)之间寻求最佳折衷，以期获得最好的分类能力。在使用前需要对支持向量机svm进行训练，本发明在图像库中获取100张云图像和100张雪图像，并提取这100张云图像的gabor特征和100张雪图像的gabor特征，从而对支持向量机svm进行训练。其中，在提取图像库中的云图像和雪图像的gabor特征时，gabor滤波器核函数设置的方向与正弦函数波长，与步骤s3中提取云雪区域对应的超像素块的gabor特征时，gabor滤波器核函数设置的方向与正弦函数波长相同。也就是说，本实施例中，支持向量机svm训练采用的特征就是上述的8个gabor核函数方向与2个正弦函数波长确定的16维gabor特征。训练时使用的核函数是高斯核函数，其余参数选用模型默认值。
需要说明的是，由于云是漂浮在空中，由汽化的水汽不断聚合而形成的，则云周围梯度过渡平滑，而雪是紧贴地表周围过渡明显，边缘纹理较大。同时汽化的云的纹理方向没有规律，而由于地面纹理的影响，雪区域纹理有明显的方向性，例如图5(a)为本实施例提供的可见光波段云区域遥感图像示意图，图5(b)为本实施例提供的可见光波段雪区域遥感图像示意图；图6(a)为图5(a)对应的梯度统计直方图；图6(b)为图5(b)对应的梯度统计直方图。基于步骤s2的获取的云雪区域，进一步提取云雪区域对应的超像素块的gabor特征。参见图7(a)、图7(b)，分别为本实施例提供的可见光波段遥感图像示意图、以及可见光波段遥感图像示意图对应的gabor特征示意图。其中，图7(a)白色部分为云雪区域，灰色部分为山区或阴影。图7(b)为gabor值图，gabor值越小，颜色越深，则图7(b)中深黑色区域为云区域，白色边界为雪区域。比对图7(a)和图7(b)，本实施例的基于空域自然场景统计和gabor特征的遥感图像云检测方法，能够准确地从可见光波段遥感图像中提取出云区域。
当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当然可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。