基于支持向量机的光学卫星遥感影像云雪雾检测方法

摘要

本发明公开了一种基于支持向量机的卫星遥感影像云、雪、雾检测方法，包括以下步骤：首先，收集不同类型的大量地物和云、雪、雾样本影像数据影像作为训练集，获得影像的灰度和纹理特征组成特征集合；通过支持向量机的方法对所有样本的特征集合进行机器学习获得云、雪、雾影像分类器。其次，使用得到的云、雪、雾影像分类器确定待测影像的类别，并进行形态学闭运算和重叠区域校正，判断遥感影像中目标区域类型；最后，重新选择训练样本获得新的影像分类器，对待测卫星遥感影像进行二次检测，并与第一次检测作比较，最终确定待测遥感影像的云、雪、雾的判定结果。实验结果表明本发明方法能够获得较高的检测精度。

说明书

技术领域

本发明属于卫星遥感影像质量检测领域，具体涉及一种基于支持向量机的卫星遥感影像云、雪、雾检测方法。

背景技术

在光学卫星遥感影像中，遥感信息常常受到云、雾天气和积雪的影响，云和积雪会覆盖影像所在区域的地表信息，雾与霾等则会使遥感影像中许多特征信息被掩盖。因此，需要对遥感影像中的云、雪、雾区域进行检测，并剔除无效信息覆盖过大的相应影像数据，从而提高光学卫星遥感影像的利用率。

目前的遥感影像云、雪、雾检测方法主要集中在对云或雾的检测、云、雾和云、雪的检测识别。主要包括阈值法，特征提取法等。遥感影像云检测的方法主要通过利用不同波段下的光谱反射率来设定光谱阈值判断是否是云，或通过提取影像特征，根据特征分类来进行云检测。遥感影像雾检测方法主要是对典型个例的研究，利用遥感数据提取特征进行大雾的监测研究。云、雪的检测方法主要利用云、雪在可见光波段特征相似而在短波红外差异较大的特点，通过构建云、雪反差增大因子识别雪，或通过对全色影像的纹理特征计算分形维数实现云、雪的识别。而云、雪、雾检测方法通常是上述方法的叠加使用。

经过对现有文献检索发现，目前的云、雪、雾检测方法有如下问题：第一、现有的方法难以同时检测云、雪、雾。检测方法受检测类型的影响，单个检测方法难以适应多种检测要求；已有的阈值法可靠性不高，由于检测结果受时空类型的影响，难以推广到较为普遍的检测上；特征提取法选取的影像特征信息不足，检测准确率不够高。第二、对云、雪、雾检测方法的检测效率不高，算法的复杂度较高，对大数据影响难以快速的检测和识别，对遥感数据源有一定要求，普适性较差。

发明内容

发明的目的在于增强遥感影像质量检查的时效性，提高遥感影像利用率，使其可以应用到资源一号、资源三号、天绘一号以及高分一号等国产卫星影像产品质量检查系统中。

为了实现这一目的，本发明提出了一种基于支持向量机的卫星遥感影像云、雪、雾检测方法，本发明技术方案的具体实现包括以下步骤：

步骤1，收集大量的云、雪、雾和地物样本影像数据；

步骤2，提取各类样本影像的灰度特征和纹理特征组成特征向量；

步骤3，利用支持向量机训练样本影像的特征向量，分别得到由决策函数构成的云影像分类器、雪影像分类器和雾影像分类器；

步骤4，对待测卫星遥感影像的原始影像进行降采样处理以获取缩略图，对缩略图进行影像切分得到子影像，计算所有子影像的灰度特征和纹理特征组成的特征向量；

步骤5，对待测卫星遥感影像的子影像进行分类，包括以下子步骤，

步骤5.1，将步骤4提取的特征向量分别输入到步骤3得到的云、雪、雾影像分类器中进行预测分类；

步骤5.2，将全部子影像按照目标区域的类型划分为云区、雾区、雪区和地物区；

步骤5.3，按照云与地物区，雾与地物区，雪与地物区划分为三幅二值化图像，其中每幅图像中的地物区取相同的零值，云、雪、雾区取不同的图像值；

步骤6，对步骤5得到的分类结果进行形态学“闭”运算；

步骤7，比较同一位置三幅二值化的图像值，获得待测卫星遥感影像中云、雪、雾的检测结果。

进一步的，所述步骤7的实现方式如下，

步骤7.1，比较同一位置三幅二值化的图像值，若三幅图像同一位置的图像值相同，则判定该位置为地物区；若同一位置存在两种相同的值，则判定该位置为第三种图像值所代表的类别区域；若同一位置三幅图像的图像值均不同，则判定该位置是存在云、雪、雾的重叠区域，并将其中存在零值的点记为二重叠区域，将不存在零值的点记为三重叠区域；

步骤7.2，重复步骤7.1，比较三幅二值化图像的所有图像值，得到云、雪、雾区和地物区以及重叠区域的判别结果，并对重叠区域进行校正；首先判断重叠区域是否包含于其它区域，若包含则将该重叠区域替换为其他区域；其次判断二重叠区域类别，若与某一确定类别区域外接，则判定该重叠区域的类别为重叠区域除去确定的外接类别后的类别，若没有则待三重叠区域类别判定完毕后确认；对于三重叠区域，若与某一确定类别区域外接，则判定该重叠区域为除去外接类别后的二重叠区域，若与二重叠区域外接，则判定该重叠区域的类别为重叠区域与其外接的二重叠区域类别不相同的类别；最后判定二重叠区域仅外接不同二重叠区域的情况，分别将这些区域判定为除去共有类别之后的类别，最终得到判定结果。

步骤7.3，将判定结果进行形态学“闭”运算，得到待测卫星遥感影像中云、雪、雾的检测结果。

进一步的，还包括步骤8，重新选取适量云与地物、雾与地物、雪与地物样本作为训练样本，重复步骤2-7对待测卫星遥感影像进行二次检测，将二次检测的检测结果与第一次检测结果进行比对，若同一位置两次检测结果相同，则确定该位置的类别为任意一次检测结果得到的类别，若同一位置两次检测结果不同，则判定该位置的类别为地物，最终得到检测结果。

进一步的，所述步骤2的实现方式如下，

步骤2.1，计算样本影像的灰度特征，包括样本影像的灰度均值、灰度方差、一阶差分和直方图信息熵；

其中灰度均值的计算公式为，

其中，f(i,j)为(i,j)处的灰度值，S＝M×N,M是样本影像的宽，N是样本影像的高；

灰度方差的计算公式为，

一阶差分的计算公式为，

直方图信息熵的计算公式为，

其中，h[g]是样本影像的直方图，h[g](i)是在某灰度级下的像素所占整个样本影像的百分比，M为最大灰度级；

步骤2.2，计算样本影像的纹理特征，包括样本影像的梯度标准差、混合熵、逆差矩、

纹理分数维；

其中梯度标准差的计算公式为，

G(i,j；d,θ)＝#{(x₁,y₁)(x₂,y₂)|f(x₁,y₁)＝i,f(x₂,y₂)＝j,|(x₁,y₁)-(x₂,y₂)|＝d,∠((x₁,y₁),(x₂,y₂))＝θ}其中d表示两个像素间的距离，θ表示像素间的方向角，f(x₁,y₁)和f(x₂,y₂)分别表示(x₁,y₁)和(x₂,y₂)的灰度值，∠表示像素点与水平位置的夹角，#表示按照集合里的限制条件得到的像素对的个数，∑#L_xL_y表示在特定的位置关系下所有像素对的个数的总和；L_g表示灰度级最大值，L表示梯度的最大值；

混合熵的计算公式为，

逆差距的计算公式为，

利用分形布朗随机场法求解样本影像的纹理分数维，则影像的分数维D的表达式为，

D＝n+1-H

其中，n指样本影像空间维数，H为自相似参数；

步骤2.3，将以上灰度特征和纹理特征组成8维特征向量。

进一步的，所述步骤3的实现方式如下，

步骤3.1，选取部分云样本和地物样本作为训练样本，将训练样本的特征向量作为训练影像分类器的训练集T＝{(x₁,y₁)...(x_i,y_i)}，i＝1…N，y_i∈ψ＝{-1，1}，其中1代表正类，代表有云区域类别，-1代表负类，代表地物区域类别，x_i∈Rⁿ，x_i为特征向量，N为样本个数；

步骤3.2，采用C-SVC模型的支持向量机构建分类超平面，计算高斯核函数，

其中，x_i与x_j分别指样本i和j的特征向量，||x_i-x_j||²为欧式距离的平方，σ为方差，i＝1…N，j＝1…N，N为样本个数；

步骤3.3，采用凸二次规划方法求解特征空间最优分类超平面的拉格朗日乘子向量，

其中，0≤α_i≤C,i＝1,2…N，为拉格朗日乘子向量，α_i和αj分别表示第i个和第j个拉格朗日乘子，x_i,x_j分别表示第i个样本和第j个样本的特征向量，y_i和y_j分别表示为第i个样本和第j个样本的类别，C为惩罚参数，N为样本个数，得到拉格朗日乘子向量的最优解为：

α^*＝(α₁^*,α₂^*,...α_i^*...α_N^*)^T

其中，α_i^*表示第i个拉格朗日乘子的最优解；

步骤3.4，求解特征空间最优分类超平面的截距，计算公式为，

其中，α_i^*为第i个拉格朗日乘子的最优解，y_i为第i个样本的类别，N为样本个数；

步骤3.5，将求得的高斯核函数，拉格朗日最优解，超平面截距代入决策函数，

步骤3.6，将剩余的云样本和地物样本作为测试样本，对决策函数进行测试，优化决策函数，同时获得相应的云影像分类器；

步骤3.7，重复步骤3.1-3.6，分别获得雪影像分类器和雾影像分类器。

进一步的，所述步骤4中，若待测遥感影像为全色影像，则直接采取降采样处理，若待测遥感影像为多光谱影像，则采取RGB三波段进行降采样。

进一步的，所述步骤6的实现方式为，选取结构形状为正方形，结构尺寸为3×3的结构元素，对三幅二值化图像分别进行膨胀运算，将得到的处理图像，再用相同的结构元素进行腐蚀运算。

与现有技术相比，本发明的优点：本发明方法可以一次训练后用于多次检测，通过大量影像训练得到影像分类器，检测时只需要再次使用即可，支持向量机算法在预测分类阶段时间复杂度低，可以快速进行区域类型的检测；经测试，本发明方法既适用于全色影像，也适用于n通道多光谱影像，利用本发明方法对资源一号02星号、资源三号、天绘一号、高分一号等多颗国产卫星遥感影像进行云检测，其准确度分别达到94.8％、96.4％、93.2％和95.2％。

附图说明

图1为本发明实施例的实现流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，此处所描述的实施示例仅用于说明和解释本发明，但并不限定本发明的保护范围。

参照图1，本发明以资源一号02C号，资源三号卫星全色影像数据以及天绘一号的多光谱遥感影像数据为例，其实现步骤如下：

步骤1，样本获取

将样本的原始影像降采样为1024×1024像素8位bmp格式缩略图，对缩略图进行影像切分。若遥感影像为全色影像，则直接采取降采样处理，若遥感影像为多光谱影像，则采取RGB三波段进行降采样。将全色影像切分为32×32的样本块，将多光谱影像切分为16×16的影像块，分别选取1500个地物样本，1000个云样本，1000个雪样本，1000个雾样本作为训练样本数据。

步骤2，特征提取：云、雪、雾和地物在灰度信息上存在许多不同之处，一般而言，全色与多光谱影像中云区域比雾区域的平均亮度要高很多，而雾区域的亮度值要大于地物区域的亮度值，云与雪区域则有相似的光谱特征，同时，云、雪、雾与地物在灰度分布以及灰度变化等方面也存在明显差异，因此利用影像的灰度特征可将目标区域影像和地物影像进行一定程度的区分。

提取样本影像的灰度特征和纹理特征矢量值形成8维特征向量，其具体实现步骤如下：

步骤2.1，计算样本影像的灰度特征

步骤2.1.1，求灰度均值，使用如下公式计算：

其中，f(i,j)为(i,j)处的灰度值，S＝M×N,M是样本影像的宽，N是样本影像的高。

步骤2.1.2，计算样本影像的灰度方差：

步骤2.1.3，计算样本影像的一阶差分：

步骤2.1.4，计算样本影像的直方图信息熵：

其中，h[g]是样本影像的直方图，h[g](i)是在某灰度级(i)下的像素所占整个影像的百分比，M为最大灰度级。

步骤2.2，计算样本影像的纹理特征：从人眼的视觉特性角度出发，卫星遥感影像中云、雪、雾的纹理信息往往比地物的纹理要单一和简单，此外影像中云、雪、雾区域的边缘也较为模糊、圆润，而地物边缘则通常较为锐利、梯度大。因此，可利用卫星遥感影像的纹理信息和梯度信息来进行云、雪、雾区域与地物信息的检测和划分。此外，云、雪、雾影像之间的纹理特征也具有明显区别。对于云样本，其纹理属于随机纹理，多变而且难以检测，表现杂乱没有规律，边缘纹理较粗且模糊；雾样本纹理则比较均匀，平滑度较好，边缘形态规则；积雪样本受地面纹理的影响，具有更好的方向性，梯度变化大。通过影像灰度和影像梯度的综合信息，对云、雪、雾表现的不同纹理特征进行分辨。

步骤2.2.1，计算出样本影像的灰度梯度共生矩阵G(i,j,d,θ)，使用如下公式：

G(i,j；d,θ)＝#{(x₁,y₁)(x₂,y₂)|f(x₁,y₁)＝i,f(x₂,y₂)＝j,|(x₁,y₁)-(x₂,y₂)|＝d,∠((x₁,y₁),(x₂,y₂))＝θ}其中d表示两个像素间的距离，θ表示像素间的方向角，(x,y)表示像素点的坐标，f(x,y)表示该点的灰度值，∠表示像素点与水平位置的夹角，#表示按照集合里的限制条件得到的像素对的个数，例如两个点像素值分别为1和2，它们之间的距离为1，θ＝0°即水平方向，在整个影像中遍历，统计符合这些条件的两个点组成的像素对有多少个。

步骤2.2.2，将灰度共生矩阵G(i,j,d,θ)归一化为H(i,j；d,θ)，其计算公式如下：

其中，∑#L_xL_y表示在特定的位置关系下(即指距离d和夹角θ相同)所有像素对的个数的总和。

步骤2.2.3，计算样本影像的梯度标准差，首先计算梯度平均：

L_g表示灰度级最大值，L表示梯度的最大值。

将梯度平均值T代入如下公式得到梯度标准差：

步骤2.2.4，计算样本影像的混合熵，计算公式如下：

步骤2.2.5，提取样本影像的局部平稳性特征，计算逆差距，其计算公式如下：

步骤2.2.6，选用分形布朗随机场法来求解样本影像的纹理分数维；

求解常数H(0<H<1)，使得分布函数F(t)满足：

F(t)为x，Δx无关的分布函数，H为自相似参数，f(x)称为关于x的实随机函数，n指样本影像空间维数，则影像的分数维D的表达式为：

D＝n+1-H

步骤3，影像分类器训练

使用支持向量机的方法来训练样本影像的特征向量，分别得到由决策函数构成的云影像分类器、雪影像分类器和雾影像分类器，其具体实现包括以下子步骤：

步骤3.1，将云样本和地物样本的80％作为训练样本，将训练样本的特征向量作为训练影像分类器的训练集T＝{(x₁,y₁)...(x_i,y_i)}，i＝1…N。其中，y_i∈ψ＝{-1，1}其中1代表正类，即有云区域类别，-1代表负类，即地物区域类别，x_i∈Rⁿ，x_i为特征向量，N为样本个数。

步骤3.2，采用C-SVC模型的支持向量机构建分类超平面，计算高斯核函数：

其中x_i与x_j分别指样本i和j的特征向量，||x_i-x_j||²为欧式距离的平方，σ为方差，σ的值一般根据实验结果恰当选取。由于在分类的过程中很难完全的将特征向量正确分类，因此σ的目的可以理解为设置一个容错范围，在范围内的忽略错误。i＝1…N，j＝1…N，N为样本数量。

步骤3.3，采用凸二次规划方法求解特征空间最优分类超平面的拉格朗日乘子向量：

其中0≤α_i≤C,i＝1,2…N，为拉格朗日乘子向量，α_i和αj分别为其中的第i个和第j个拉格朗日乘子，x_i,x_j分别表示第i个样本和第j个样本的特征向量，y_i和y_j分别表示为第i个样本和第j个样本的类别，C为惩罚参数，N为样本个数，求解出拉格朗日乘子向量的最优解为：

α^*＝(α₁^*,α₂^*,...α_i^*...α_N^*)^T

其中α_i^*表示第i个拉格朗日乘子的最优解。

步骤3.4，求解特征空间最优分类超平面的截距，计算公式：

其中，α_i^*为第i个拉格朗日乘子的最优解，y_i为第i个样本的正,负类别，N为样本个数。

步骤3.5，将前述求得的高斯核函数，拉格朗日最优解，超平面截距代入决策函数：

步骤3.6，将云样本和地物样本的20％作为测试样本集，对决策函数进行测试，优化决策函数，获得相应的云影像分类器。

步骤3.7，重复步骤3.1-3.6，分别获得雪影像分类器和雾影像分类器。

步骤4，待测影像特征提取

将待测的原始影像降采样为1024×1024像素8位bmp格式缩略图，若遥感影像为全色影像，则直接采取降采样处理，若遥感影像为多光谱影像，则采取RGB三波段进行降采样，对缩略图进行影像切分以获取其1024幅32×32像素子影像，提取所有子影像的特征向量，包括灰度特征和纹理特征矢量，具体提取可通过步骤2实现。

步骤5，待测影像分类

步骤5.1，将步骤4中提取的特征向量分别输入到步骤3得到的对应的云、雪、雾影像分类器中进行预测分类，通过决策函数对特征向量进行分类。

步骤5.2，重复执行步骤3直到所有子影像都被分类，将全部子影像按照目标区域的类别划分为云区、雾区、雪区和非云雪雾区域(即地物区)。

步骤5.3，按照云与地物区，雾与地物区，雪与地物区划分为三幅二值化图像，其中各幅图像之间的地物区取相同的零值，而云、雪、雾区域之间则取不同的图像值。

步骤6，形态学闭运算

选取结构形状为正方形，结构尺寸为3×3的结构元素，对三幅二值化图像分别进行膨胀运算，将得到的处理图像，再用相同的结构元素进行腐蚀运算，将云、雪、雾区域连成一片，最终消除边缘的噪声区域。

步骤7，重叠区域校正

步骤7.1，比较同一位置三幅二值化的图像值：若三幅图像同一位置的图像值相同，则判定该位置为地物区；若存在两种相同的值，则表明该位置为第三种图像值所代表的类别区域；若图像值不同，则表明该子影像存在云、雪、雾的重叠区域，将其中存在零值的点记为二重叠区域，将不存在零值的点记为三重叠区域。

步骤7.2，重复步骤7.1，比较三幅二值化图像的所有图像值，得到云、雪、雾区域和地物区以及重叠区域的判别结果,并对重叠区域进行校正；首先判断重叠区域是否包含于其它区域，若包含则将该重叠区域替换为其他区域(包括重叠区域和确定类别区域)；其次判断二重叠区域类别，若与某一确定类别区域外接，则判定该重叠区域的类别为重叠区域除去确定的外接类别后的类别，若没有则待三重叠区域类别判定完毕后确认；对于三重叠区域，若与某一确定类别区域外接，则判定该重叠区域为除去外接类别后的二重叠区域，若与二重叠区域外接，则判定该重叠区域的类别为重叠区域与其外接的二重叠区域类别不相同的类别；最后判定二重叠区域仅外接不同二重叠区域的情况，分别将这些区域判定为除去共有类别之后的类别，最终得到判定结果。例如云雪重叠区域外围如果是确定的雾区域，则该云雪区域判定为雾区域；如果云雪区域被云雾区域包围，则判定该云雪区域为云雾区域；如果云雪区域与确定的云区域外接，则判断该区域为雪区域；如果云雪雾区域确定的云域外接，则判定该区域为雪雾区域；如果云雪雾区域与云雾区域外接，则判断该区域为确定的雪区域；如果云雪区域与云雾区域外接，则判断该云雪区域与云雾区域分别为雪区域和雾区域。

步骤7.3，将判定结果进行形态学“闭”运算，其运算方法如步骤6所述，得到最终的云、雪、雾检测结果。

步骤8，二次检测

重新选取各500个云与地物、雾与地物、雪与地物样本制作支持向量机分类器，地物选取高亮样本。对待测样本进行“二次检测”，将二次检测的检测结果与第一次检测结果进行比对，若同一位置两次检测结果相同，则确定该位置的类别为任意一次检测结果得到的类别，若同一位置两次检测结果不同，则判定该位置的类别为地物，最终得到检测结果。例如，如果二检的结果是云或雪，一检的结果是雾，该区域都判定为地物，只有一检和二检的结果都是云，才能确定该区域为云。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。