CN105354865A - 多光谱遥感卫星影像自动云检测方法及系统 - Google Patents

Abstract

本发明提供一种多光谱遥感卫星影像自动云检测方法及系统，包括数据准备，云层粗提取，在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，将原影像的强度信息作为向导图，在误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取。本发明技术方案真正实现了云的自动、快速、准确检测。<b />

Description

多光谱遥感卫星影像自动云检测方法及系统

技术领域

本发明属于测绘科学与技术领域，涉及一种基于阈值分割与纹理信息的多光谱遥感卫星影像自动云检测技术。

背景技术

卫星影像成像过程中，由于受到大气影响，会存在云遮挡，导致影像质量受到严重影响。厚云对地面光谱信息造成无法弥补的破坏，薄云虽未完全遮蔽地物信息，但仍会造成原地物光谱失真，严重影响地物判读、测绘卫星产品生产等后续工作的进行。因此，实现生产前准确云检测具有重要的实际生产意义。

目前主要利用云的光谱、频率、纹理等特性，结合阈值法、支持向量机法、聚类法等进行检测。众多算法中，光谱结合阈值的方法最为普遍，主要利用云在可见光、近红外等波段具有强反射的特性，多以像素为单位，具有数据质量要求低、计算速度快、适用性强等优点，但是存在阈值敏感问题，同一个卫星数据因时间、天气、日照、湿度等原因可能检测阈值将发生巨大变化，加大了光谱结合阈值的云检测方法局限性。频率结合阈值的方法主要利用云的低频率特性，通过小波分析、傅里叶变换等方法获取影像低频数据来进行云检测，但受到低频信息的干扰，通常采用多层小波变换排除干扰，这大大降低了云检测效率。纹理特征法主要利用云与地面纹理特征差异，一般来说云纹理细节明显少于地面，常以分块子图为单位，结合分形维数和角二阶矩法进行纹理特征计算，但是该类方法需要足够的训练样本以获得可靠的云特征区间才能保证分类的精度，导致效率低下。

发明内容

本发明的目的是提供一种快速、准确、全自动化的基于阈值分割与纹理信息的多光谱遥感卫星影像自动云检测，它能够克服上述现有云检测方法技术的不足，满足高分辨率遥感卫星影像云检测应用的需求。

为达到上述目的，本发明提供的技术方案提供一种多光谱遥感卫星影像自动云检测方法，包括以下步骤：

步骤1，数据准备，包括获取待进行云检测的卫星影像，卫星影像中包含红、绿、蓝、近红外四个波段的数据；

步骤2，云层粗提取，包括将原影像变换到HIS色彩空间，将强度信息和饱和度信息分别线性拉伸至[01]，然后计算得到基底影像，判断基底像素值大于基底阈值且近红外像素值大于近红外阈值且色调像素值小于色调阈值的像元为粗检云对象，否则为非云对象，得到云层粗检结果；

将强度信息和饱和度信息线性拉伸至[01]后，利用下述公式计算得到基底影像，

其中，I’代表原影像中对应像素的强度值线性拉伸后结果，S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，J代表生成的基底影像对应像素值；

步骤3，云层误差剔除，包括在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，

步骤4，云层精提取，包括将步骤2变换所得原影像的强度信息作为向导图，在步骤3进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，包括以下子步骤，

步骤4.1，初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为k1，设定新增面积阈值，迭代次数阈值D；

步骤4.2，设在步骤3进行误差剔除后粗检结果为影像I_R，在影像I_R中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为种子点；

步骤4.3，若满足m>T且d<D则进入步骤4.4，若满足m<T且k＝k1则进入步骤4.5，否则将影像I_R中面积小于预设面积阈值的云区域清除，然后结束流程；

步骤4.4，在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，返回步骤4.2；

步骤4.5，在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，返回步骤4.2；

其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。

而且，步骤2中，所述基底阈值采用以下方式确定，

首先，计算Otsu阈值采用如下公式，

T＝MAX{U_t0×(V_to-V_tf)²+U_t1×(V_t1-V_tf)²}

其中，U_t0表示当阈值为t时，前景像元占总像元数比例，V_t0表示前景像元平均灰度，U_t1表示背景像元占总像元数比例，V_t1表示背景像元平均灰度，V_tf表示影像总平均灰度，V_tf＝U_t0×V_t0+U_t1×V_t1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t作为Otsu阈值T_J；

设云检测基底灰度阈值变化空间为[T_aT_b]，T_a、T_b分别为云检测基底灰度阈值变化下限和上限，按下式得到修正后的Otsu阈值作为基底阈值，

$T_{JF}=\left\{\begin{array}{ccc}{T}_a,& if& T_J<T_a\\ T_J,& if& T_a\&le;T_J\&le;T_b\\ T_b,& if& T_J>T_b\end{array}\right.$

其中T_J为Otsu阈值，T_JF为修正后的Otsu阈值。

而且，步骤3中，

所述直方图均衡化的计算公式如下，

$\begin{array}{cc}{S}_k=\frac{\left({\mathrm{\&Sigma;}}_{j=0}^k{n}_j\right)}{N}\&times;L& k\&Element;\left(\begin{array}{cc}0& L\end{array}\right)\end{array}$

公式中，S_k表示均衡化后的灰度值，表示灰度色阶0～k像素个数，k属于0到L；n_j表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶；

将直方图均衡化之后的强度灰度值进行双边滤波处理，包括将直方图均衡化之后强度S_k按下式得到滤波影像I’，

$I^{,}(i,j)=\frac{{\&Sigma;}_{\left(i,j\right)\&Element;{\mathrm{\&Omega;}}_{x,y}}w\left(i,j\right)S_k\left(i,j\right)}{{\&Sigma;}_{\left(i,j\right)\&Element;{\mathrm{\&Omega;}}_{x,y}}w\left(i,j\right)}$

w(i,j)＝W_s(i,j)×W_r(i,j)

$W_s(i,j)=e^{\frac{-\left|\right|(i,j)-(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_s}^2}}$

$W_r(i,j)=e^{\frac{-\left|\right|s_k(i,j)-s_k(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_r}^2}}$

D＝|S_k-I’|

其中，Ω_x,y表示中心点(x，y)的M×M大小的领域，对该领域内每一个像素点(i,j)均衡化后强度值S_k(i,j)，记w(i,j)为计算权重，W_s(i,j)为空间距离因子，W_r(i,j)为亮度相似度因子，σ_s和σ_r表示分别相应的高斯函数标准偏差，表示均衡化影像中强度最大值，e表示(i,j)为待滤波中心像素坐标，(i’,j’)为M*M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。

本发明还相应提供一种多光谱遥感卫星影像自动云检测系统，包括以下模块：

数据准备模块，用于获取待进行云检测的卫星影像，卫星影像中包含红、绿、蓝、近红外四个波段的数据；

云层粗提取模块，用于将原影像变换到HIS色彩空间，将强度信息和饱和度信息分别线性拉伸至[01]，然后计算得到基底影像，判断基底像素值大于基底阈值且近红外像素值大于近红外阈值且色调像素值小于色调阈值的像元为粗检云对象，否则为非云对象，得到云层粗检结果；

将强度信息和饱和度信息线性拉伸至[01]后，利用下述公式计算得到基底影像，

其中，I’代表原影像中对应像素的强度值线性拉伸后结果，S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，J代表生成的基底影像对应像素值；

云层误差剔除模块，用于在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，

云层精提取模块，用于将云层粗提取模块变换所得原影像的强度信息作为向导图，在云层误差剔除模块进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，包括以下子模块，

第一子模块，用于初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为k1，设定新增面积阈值，迭代次数阈值D；

第二子模块，用于设在云层误差剔除模块进行误差剔除后粗检结果为影像I_R，在影像I_R中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为种子点；

第三子模块，用于若满足m>T且d<D则命令第四子模块工作，若满足m<T且k＝k1则命令第五子模块工作，否则将影像I_R中面积小于预设面积阈值的云区域清除，然后结束公子；

第四子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，命令第二子模块工作；

第五子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，命令第二子模块工作；

其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。

而且，云层粗提取模块中，所述基底阈值采用以下方式确定，

首先，计算Otsu阈值采用如下公式，

T＝MAX{U_t0×(V_to-V_tf)²+U_t1×(V_t1-V_tf)²}

其中，U_t0表示当阈值为t时，前景像元占总像元数比例，V_t0表示前景像元平均灰度，U_t1表示背景像元占总像元数比例，V_t1表示背景像元平均灰度，V_tf表示影像总平均灰度，V_tf＝U_t0×V_t0+U_t1×V_t1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t作为Otsu阈值T_J；

设云检测基底灰度阈值变化空间为[T_aT_b]，T_a、T_b分别为云检测基底灰度阈值变化下限和上限，按下式得到修正后的Otsu阈值作为基底阈值，

$T_{JF}=\left\{\begin{array}{ccc}{T}_a,& if& T_J<T_a\\ T_J,& if& T_a\&le;T_J\&le;T_b\\ T_b,& if& T_J>T_b\end{array}\right.$

其中T_J为Otsu阈值，T_JF为修正后的Otsu阈值。

而且，云层误差剔除模块中，

所述直方图均衡化的计算公式如下，

$\begin{array}{cc}{S}_k=\frac{\left({\mathrm{\&Sigma;}}_{j=0}^k{n}_j\right)}{N}\&times;L& k\&Element;\left(\begin{array}{cc}0& L\end{array}\right)\end{array}$

公式中，S_k表示均衡化后的灰度值，表示灰度色阶0～k像素个数，k属于0到L；n_j表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶；

将直方图均衡化之后的强度灰度值进行双边滤波处理，包括将直方图均衡化之后强度S_k按下式得到滤波影像I’，

$I^{,}(i,j)=\frac{{\&Sigma;}_{\left(i,j\right)\&Element;{\mathrm{\&Omega;}}_{x,y}}w\left(i,j\right)S_k\left(i,j\right)}{{\&Sigma;}_{\left(i,j\right)\&Element;{\mathrm{\&Omega;}}_{x,y}}w\left(i,j\right)}$

w(i,j)＝W_s(i,j)×W_r(i,j)

$W_s(i,j)=e^{\frac{-\left|\right|(i,j)-(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_s}^2}}$

$W_r(i,j)=e^{\frac{-\left|\right|s_k(i,j)-s_k(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_r}^2}}$

D＝|S_k-I’|

其中，Ω_x,y表示中心点(x，y)的M×M大小的领域，对该领域内每一个像素点(i,j)均衡化后强度值S_k(i,j)，记w(i,j)为计算权重，W_s(i,j)为空间距离因子，W_r(i,j)为亮度相似度因子，σ_s和σ_r表示分别相应的高斯函数标准偏差，表示均衡化影像中强度最大值，e表示(i,j)为待滤波中心像素坐标，(i’,j’)为M*M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。

本发明所采用多光谱阈值与纹理信息相结合的方法进行云检测，首先采用改进的颜色模型转换方法，将影像由RGB转换至HIS颜色空间。利用影像强度信息与饱和度信息结合Otsu阈值分割生成基底图，并结合影像近红外信息与色调信息对基底图进行优化，生成修正图。然后利用直方图均衡化与双边滤波结合二维Otsu分割方法提取地物的纹理信息，对修正图进行粗差剔除，生成云种子图。最后将强度信息作为向导结合种子图信息进行云精确提取。本发明主要应用于卫星影像云检测，保证卫星影像产品生产等后续工作的顺利进行、大范围高质量无云遥感影像合成等。本发明技术方案具有以下优点：

1)利用改进的HIS模型，将原影像变换到HIS色彩空间，充分利用了云的光谱特性，大大降低了薄云的检测难度。

2)利用直方图均衡化和双边滤波相结合的方法进行纹理信息提取，充分体现了云层粗纹理特征，能够将云与高相似度下垫面有效区分。

3)二维Otsu阈值分割法，提高了阈值的适用性，并提高了纹理检测精度。

4)边缘种子膨胀，大大提高了云层检测精度。

附图说明

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

图2是本发明实施例的云层基底像素值统计示意图，其中图2(a)为灰度直方图，图2(b)为经验累积分布图。

图3是本发明实施例的二维Otsu阈值生成原理示意图，其中图3(a)为一维Otsu阈值示意图，图3(b)为二维Otsu阈值示意图。

图4是本发明实施例的膨胀操作原理示意图。

图5是本发明实施例的云层精提取流程图。

具体实施方式

以下结合实施例和附图详细说明本发明技术方案。

本发明提供一种快速、准确、全自动化的基于阈值分割与纹理信息的多光谱遥感卫星影像自动云检测方法，主要利用影像的光谱以及纹理信息进行云提取，首先利用光谱信息进行云结果初步检测，然后利用纹理信息对初检结果中的误差进行剔除，最后回头利用光谱信息对云层进行精确检测。输入数据为能提供红、绿、蓝、近红外波段数据的各类高分辨遥感卫星影像。参见图1，实施例所提供具体实现方法包含以下步骤：

步骤1：数据准备。获取待进行云检测的卫星影像，该影像中必须包含红、绿、蓝、近红外四个波段的数据。

本方法适用于所有能够提供红、绿、蓝、近红外四个波段的高分辨率遥感卫星影像数据。研究发现，RGB彩色模型对厚云检测效果较好，但难以识别薄云，并且人的视觉对亮度的敏感程度远强于对颜色浓淡的敏感程度，因此HIS色彩空间比RGB色彩空间更符合人的视觉特性。因此，后续在步骤2通过改进的HIS模型仿真人眼感知色彩的特性，充分利用云的光谱特性，大大降低了薄云的检测难度。

步骤2：云层粗提取。将原影像变换到HIS色彩空间，将强度信息I、饱和度信息S线性拉伸至[01]。然后计算得到基底影像。具体实施时，近红外阈值、色调阈值可根据统计结果预设，例如从高分一号、资源三号等卫星影像中，选取500个云对象(包含薄云、厚云、卷云等)，统计发现，云近红外值均大于350，色调值均小于120。通过带限定条件的Otsu阈值(作为基底阈值)结合近红外、色调信息的限定条件对基底影像进行分割，即基底像素值大于基底阈值(J大于T_JF)并且近红外像素值大于近红外阈值并且色调像素值小于色调阈值的像元即为粗检云对象，而不符合任意某一条件的像元，即为非云对象，以此得到云层粗检结果，即对基底图去除无云区域的修正图。

本发明涉及的RGB转HIS色彩通道计算公式，如以下所示：

$\left\{\begin{array}{c}{R}^{,}=min(R,G,B)\\ G^{,}=mid(R,G,B)\\ B^{,}=\max (R,G,B)\end{array}\right.$

I＝1/3(R’+G’+B’)

$S=1-\&lsqb;\frac{1}{3(R^{,}+G^{,}+B^{,})}\&rsqb;\&times;min(R^{,},G^{,},B^{,})$

$H^0=arccos\{\sqrt{2}/4\&times;\&lsqb;(2R^{,}-\sqrt{3}{G}^{,})+(\sqrt{2}{R}^{,}-B^{,})\&rsqb;/{\&lsqb;{(\sqrt{2}{R}^{,}-\sqrt{6}{G}^{,}/2)}^2+{(\sqrt{2}{R}^{,}-B^{,})}^2\&rsqb;}^{0.5}\}$

$H=\left\{\begin{array}{cc}0,& ifS=0\\ 360-H^0\&times;360/2\mathrm{\&pi;},& if(\sqrt{2}{B}^{,}>\sqrt{3}{G}^{,})\\ H^0\&times;360/2\mathrm{\&pi;},& otherwise\end{array}\right.$

式中R、G、B为原始影像色彩灰度值，R’为R、G、B中最小值，G’为R、G、B中间值，B’为R、G、B中最大值，min(R’,G’,B’)为R’,G’,B’中最小值，H⁰为影像色调值的初始值，H为影像色调值，S为影像饱和度信息，I为影像强度信息(即亮度)。

本发明使用影像强度、饱和度、色调、近红外信息进行云层粗检测，其策略如下：

将原影像变换到HIS色彩空间，并将I、S线性拉伸至[01](即变换到0-1之间)，并利用下述公式计算得到基底影像。

$J=\frac{I^{,}+\mathrm{\&tau;}}{S^{,}+\mathrm{\&tau;}}$

在基于像素进行计算时，公式中影像强度信息I’代表原影像中对应像素的强度值线性拉伸后结果，影像饱和度信息S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，具体实施时可以由本领域技术人员预设为大于0的值，实施例中取值为1.0，J代表生成的基底影像对应像素值，大小属于0-255。

利用带限定条件的Otsu阈值法，对基底影像分割，生成云层粗检测结果。计算公式如下：

T＝MAX{U_t0×(V_to-V_tf)²+U_t1×(V_t1-V_tf)²}

其中U_t0表示当阈值为t时，前景像元占总像元数比例，V_t0表示前景像元平均灰度，U_t1表示背景像元占总像元数比例，V_t1表示背景像元平均灰度，V_tf表示影像总平均灰度，V_tf＝U_t0×V_t0+U_t1×V_t1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t即为Otsu阈值T_J。Otsu阈值大体可将云层与下垫面分割开，但往往有分割阈值过大，薄云漏检；分割阈值过小误检率过高等现象存在，通过带约束的Otsu阈值能有效解决这一问题。

根据云检测基底灰度阈值变化空间，可以对阈值进行修正，设云检测基底灰度阈值变化空间为[T_aT_b]，T_a、T_b分别为云检测基底灰度阈值变化下限和上限，按下式得到修正阈值，

$T_{JF}=\left\{\begin{array}{ccc}{T}_a,& if& T_J<T_a\\ T_J,& if& T_a\&le;T_J\&le;T_b\\ T_b,& if& T_J>T_b\end{array}\right.$

具体实施时，云检测基底灰度阈值变化下限和上限可根据统计结果预设，例如从高分一号、资源三号等卫星影像中，选取500个云对象(包含薄云、厚云、卷云等)，统计其灰度直方图和经验累积分布图如图2(a)、图2(b)所示。

根据图2分析，云层基底像素值大体分布在[80255]之间，累积频率在130开始呈现快速上升的趋势，表征当某基底像素灰度值为130时，为云的概率极大，此时如果计算Otsu阈值大于130，那么漏掉薄云的可能性将随之迅速增大。因此，实验中设定云检测基底灰度阈值变化空间为[80130]，即按下式得到修正阈值。

$T_{JF}=\left\{\begin{array}{ccc}80,& if& T_J<80\\ T_J,& if& 80\&le;T_J\&le;130\\ 130,& if& T_J>130\end{array}\right.$

其中T_J为Otsu阈值，T_JF为修正后的Otsu阈值。

步骤3：云层误差剔除。在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方法进行纹理信息提取。然后利用二维Otsu阈值T_d(即在一维Otsu阈值基础上，将满足阈值条件的像素作为下一轮直方图统计数据来源，再次计算Otsu阈值)对纹理信息图进行分割，利用分割后的细节图对步骤2所得粗检结果进行误差剔除，即根据细节图进行粗差剔除得到云种子图。

本发明首先使用直方图均衡化突出隐含有纹理细节的图像。通过控制图像灰度级概率密度函数，改变图像的灰度层次。均衡化之后，山区灰度跳变幅度变大，城区呈现明显颗粒感，突出了建筑物等细节。而云层只是对比度增强，纹理却没有变细，灰度空间分布仍保持较好的连续性。

直方图均衡化计算公式如下：

$\begin{array}{cc}{S}_k=\frac{\left({\mathrm{\&Sigma;}}_{j=0}^k{n}_j\right)}{N}\&times;L& k\&Element;\left(\begin{array}{cc}0& L\end{array}\right)\end{array}$

公式中，S_k表示均衡化后的灰度值，表示灰度色阶0～k像素个数，k属于0到L；n_j表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶。

然后将直方图均衡化之后的强度灰度值进行双边滤波处理，双边滤波对像素同时使用空间和边界滤波，在对图像进行光滑的同时能保持图像边缘特征。对于细节信息较少的卫星影像，通常需要多次选权双边滤波提取地物纹理，但对计算机内存及计算速度是一个极大考验。文中直方图均衡化处理，有效强化了地物纹理，仅需一次双边滤波就能有效提取地物纹理。将直方图均衡化之后强度S_k按下式得到滤波影像I’：

$I^{,}(i,j)=\frac{{\&Sigma;}_{\left(i,j\right)\&Element;{\mathrm{\&Omega;}}_{x,y}}w\left(i,j\right)S_k\left(i,j\right)}{{\&Sigma;}_{\left(i,j\right)\&Element;{\mathrm{\&Omega;}}_{x,y}}w\left(i,j\right)}$

w(i,j)＝W_s(i,j)×W_r(i,j)

$W_s(i,j)=e^{\frac{-\left|\right|(i,j)-(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_s}^2}}$

$W_r(i,j)=e^{\frac{-\left|\right|s_k(i,j)-s_k(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_r}^2}}$

D＝|S_k-I’|

式中Ω_x,y表示中心点(x，y)的M×M(M为奇数)大小的领域。对该领域内每一个像素点(i,j)均衡化后强度值S_k(i,j)，记w(i,j)为计算权重，W_s(i,j)为空间距离因子，W_r(i,j)为亮度相似度因子，σ_s和σ_r表示分别相应的高斯函数标准偏差，具体实施时，本领域技术人员可自行预设取值，实施例中σ_s取值为2，σ_r取值为 表示均衡化影像中强度最大值，e表示数学常量2.718281，是自然对数函数的底数，(i,j)为待滤波中心像素坐标，(i’,j’)为M×M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。

采用二维Otsu阈值对细节图进行二值分割，即在一维Otsu阈值基础上，将满足阈值条件的像素作为下一轮直方图统计数据来源，代入Otsu阈值计算公式，再次计算Otsu阈值，过程如图3所示，如图3(a)中一维Otsu阈值Otsu_1＝22，如图3(b)中二维Otsu阈值Otsu_2＝7。具体实施时，可按照获取T_J方式执行一遍，将高于阈值的像素剔除，低于阈值条件的像素保留。根据保留下来的像素及相应统计直方图数据代入Otsu阈值计算公式，再次计算T_J。例如U_t0表示当阈值为t时前景像元占总像元数比例，再次计算时采用保留的像素相应统计直方图中所有大于t灰度的直方图部分的比例。

二维Otsu阈值分割在一定程度上提高阈值的适用性，并提高了检测精度。随后，为剔除弱噪声，对分割结果进行一次3*3窗口的膨胀操作获取二值细节图。用二值细节图对粗检结果进行粗差剔除，计算公式如下：

I_R＝I_X∩I_D

式中I_R表示粗差剔除后粗检结果影像，即为云种子图，I_X表示粗检结果影像，I_D表示二值细节图。

其中，3*3窗口的膨胀操作原理如图4所示，将大小为3*3的模板与原始图卷积，若判断像素与其八领域均为检测目标则保留，否则剔除。

步骤4：云层精提取。将步骤2所得原始影像强度信息作为向导图，在步骤3进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，得到最终的云检结果。当k＝0.8‰时，能够将漏检厚云提取出来。而先将厚云向周围薄云过渡，调整灰度阈值后，能够有效将厚云周边薄云提取出来。

将步骤2所得RGB转HIS所得影像强度信息I作为向导图，在粗差剔除后粗检结果影像I_R上利用边缘种子膨胀的方式进行精确提取，实施例的云层精确提取流程如图5所示，步骤如下：

步骤4.1：初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为0.8‰，设定新增面积阈值T＝200，迭代次数阈值D＝10；具体实施时，k、T、D可由用户自行预设，例如实施例采用的实验统计最佳设定值；

步骤4.2：粗差剔除后粗检结果影像I_R中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为待处理的种子点；

步骤4.3：若满足条件一：m>T且d<D则进入步骤4.4；若满足条件二：m<T且k＝0.8‰则进入步骤4.5；若均不满足则将影像I_R中面积较小云区域清除，具体实施时，本领域技术人员可以自行预设面积阈值，例如将影像I_R中面积小于120个像素的云区域，然后结束流程；

步骤4.4：在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m、d(d＝d+1)值，返回步骤4.2；

步骤4.5：在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k为10％作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值。所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k为1.2‰，返回步骤4.2。

在步骤4.5中，通过调整k值，首先k取极小值k1＝0.8‰，将厚云周边的漏检的同类厚云搜索出来；然后将k调整为极大值k2＝10％，等同于强行将检测结果向外扩充一圈，将厚云向周边的漏检薄云过度；最后将k调整为较小值k3＝1.2‰，再次搜索，可将薄云周边漏检的同类薄云搜索出来。具体实施时，k1、k2、k3的取值可由本领域技术人员预设为经验值，建议k1,k3一般不超过1.5％,k2不超过15％。

从理论上分析，该方法以光谱结合纹理的方法进行云检测。但不同于传统阈值法的是：其一：传统光谱结合阈值法，多采用RGB彩色模型，而本方法通过改进的HIS模型，能够有效降低薄云的检测难度；其二：该方法不拘泥于阈值的选取，通过带限定条件的Otsu阈值可保证厚、薄云被提取的同时，减少同类地物的干扰；其三：阈值分割过程中，引入近红外波段信息，在区分云与高亮的河流、积雪等地物时取得理想的效果。亦不同于传统纹理分析法的是：采用直方图均衡化与双边滤波相结合提取影像纹理信息的方法，有效的弥补了传统分形维数和角二阶矩计算精度低、效率低下的不足，且对厚云、薄云均能取得理想效果。

利用本发明进行试验，包括全自动的对资源三号、高分一号、高分二号等高分辨率卫星影像进行高效率云检测，可验证分析该方法有以下优点：

无需任何先验知识，全自动化处理；

算法复杂度适中，计算量小，计算速度快；

检测精度高，可在高精确度检测出厚云、薄云的前提下，有效排除高亮度水域、房屋、裸地等干扰；

算法适用性广，可用于能提供红、绿、蓝、近红外波段数据的各类卫星影像的云检测。

本发明技术方案可采用计算机软件方式支持自动运行流程，具体实施时也可采用模块化方式提供相应系统。本发明实施例提供一种多光谱遥感卫星影像自动云检测系统，包括以下模块：

数据准备模块，用于获取待进行云检测的卫星影像，卫星影像中包含红、绿、蓝、近红外四个波段的数据；

云层粗提取模块，用于将原影像变换到HIS色彩空间，将强度信息和饱和度信息分别线性拉伸至[01]，然后计算得到基底影像，判断基底像素值大于基底阈值且近红外像素值大于近红外阈值且色调像素值小于色调阈值的像元为粗检云对象，否则为非云对象，得到云层粗检结果；

将强度信息和饱和度信息线性拉伸至[01]后，利用下述公式计算得到基底影像，

其中，I’代表原影像中对应像素的强度值线性拉伸后结果，S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，J代表生成的基底影像对应像素值；

云层误差剔除模块，用于在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，

云层精提取模块，用于将云层粗提取模块变换所得原影像的强度信息作为向导图，在云层误差剔除模块进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，包括以下子模块，

第一子模块，用于初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为k1，设定新增面积阈值，迭代次数阈值D；

第二子模块，用于设在云层误差剔除模块进行误差剔除后粗检结果为影像I_R，在影像I_R中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为种子点；

第三子模块，用于若满足m>T且d<D则命令第四子模块工作，若满足m<T且k＝k1则命令第五子模块工作，否则将影像I_R中面积小于预设面积阈值的云区域清除，然后结束公子；

第四子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，命令第二子模块工作；

第五子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，命令第二子模块工作；

其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。

各模块具体实现可参加相应步骤，本发明不予赘述。

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

Claims (6)

1.一种多光谱遥感卫星影像自动云检测方法，其特征在于，包括以下步骤：

步骤1，数据准备，包括获取待进行云检测的卫星影像，卫星影像中包含红、绿、蓝、近红外四个波段的数据；

步骤2，云层粗提取，包括将原影像变换到HIS色彩空间，将强度信息和饱和度信息分别线性拉伸至[01]，然后计算得到基底影像，判断基底像素值大于基底阈值且近红外像素值大于近红外阈值且色调像素值小于色调阈值的像元为粗检云对象，否则为非云对象，得到云层粗检结果；

将强度信息和饱和度信息线性拉伸至[01]后，利用下述公式计算得到基底影像，

其中，I’代表原影像中对应像素的强度值线性拉伸后结果，S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，J代表生成的基底影像对应像素值；

步骤3，云层误差剔除，包括在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，

步骤4，云层精提取，包括将步骤2变换所得原影像的强度信息作为向导图，在步骤3进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，包括以下子步骤，

步骤4.1，初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为k1，设定新增面积阈值，迭代次数阈值D；

步骤4.2，设在步骤3进行误差剔除后粗检结果为影像I_R，在影像I_R中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为种子点；

步骤4.3，若满足m>T且d<D则进入步骤4.4，若满足m<T且k＝k1则进入步骤4.5，否则将影像I_R中面积小于预设面积阈值的云区域清除，然后结束流程；

步骤4.4，在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，返回步骤4.2；

步骤4.5，在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，返回步骤4.2；

其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。

2.根据权利要求1所述多光谱遥感卫星影像自动云检测方法，其特征在于：步骤2中，所述基底阈值采用以下方式确定，

首先，计算Otsu阈值采用如下公式，

T＝MAX{U_t0×(V_to-V_tf)²+U_t1×(V_t1-V_tf)²}

其中，U_t0表示当阈值为t时，前景像元占总像元数比例，V_t0表示前景像元平均灰度，U_t1表示背景像元占总像元数比例，V_t1表示背景像元平均灰度，V_tf表示影像总平均灰度，V_tf＝U_t0×V_t0+U_t1×V_t1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t作为Otsu阈值T_J；

设云检测基底灰度阈值变化空间为[T_aT_b]，T_a、T_b分别为云检测基底灰度阈值变化下限和上限，按下式得到修正后的Otsu阈值作为基底阈值，

$T_{JF}=\left\{\begin{array}{cc}{T}_a,& if{T}_J<T_a\\ T_J,& if{T}_a\&le;T_J\&le;T_b\\ T_b,& if{T}_J>T_b\end{array}\right.$

其中T_J为Otsu阈值，T_JF为修正后的Otsu阈值。

3.根据权利要求1或2所述多光谱遥感卫星影像自动云检测方法，其特征在于：步骤3中，

所述直方图均衡化的计算公式如下，

$S_k=\frac{\left({\mathrm{\&Sigma;}}_{j=0}^k{n}_j\right)}{N}\&times;L,k\&Element;\left(\begin{array}{cc}0& L\end{array}\right)$

公式中，S_k表示均衡化后的灰度值，表示灰度色阶0～k像素个数，k属于0到L；n_j表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶；

将直方图均衡化之后的强度灰度值进行双边滤波处理，包括将直方图均衡化之后强度S_k按下式得到滤波影像I’，

$I^{,}(i,j)=\frac{{\mathrm{\&Sigma;}}_{(i,j)\&Element;{\mathrm{\&Omega;}}_{x,y}}w(i,j)S_k(i,j)}{{\mathrm{\&Sigma;}}_{(i,j)\&Element;{\mathrm{\&Omega;}}_{x,y}}w(i,j)}$

w(i,j)＝W_s(i,j)×W_r(i,j)

$W_S(i,j)=e^{\frac{-\left|\right|(i,j)-(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_S}^2}}$

$W_r(i,j)=e^{\frac{-\left|\right|S_k(i,j)-S_k(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_r}^2}}$

D＝|S_k-I’|

其中，Ω_x,y表示中心点(x，y)的M×M大小的领域，对该领域内每一个像素点(i,j)均衡化后强度值S_k(i,j)，记w(i,j)为计算权重，W_s(i,j)为空间距离因子，W_r(i,j)为亮度相似度因子，σ_s和σ_r表示分别相应的高斯函数标准偏差，表示均衡化影像中强度最大值，e表示(i,j)为待滤波中心像素坐标，(i’,j’)为M*M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。

4.一种多光谱遥感卫星影像自动云检测系统，其特征在于，包括以下模块：

数据准备模块，用于获取待进行云检测的卫星影像，卫星影像中包含红、绿、蓝、近红外四个波段的数据；

云层粗提取模块，用于将原影像变换到HIS色彩空间，将强度信息和饱和度信息分别线性拉伸至[01]，然后计算得到基底影像，判断基底像素值大于基底阈值且近红外像素值大于近红外阈值且色调像素值小于色调阈值的像元为粗检云对象，否则为非云对象，得到云层粗检结果；

将强度信息和饱和度信息线性拉伸至[01]后，利用下述公式计算得到基底影像，

其中，I’代表原影像中对应像素的强度值线性拉伸后结果，S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，J代表生成的基底影像对应像素值；

云层误差剔除模块，用于在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，

云层精提取模块，用于将云层粗提取模块变换所得原影像的强度信息作为向导图，在云层误差剔除模块进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，包括以下子模块，

第一子模块，用于初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为k1，设定新增面积阈值，迭代次数阈值D；

第二子模块，用于设在云层误差剔除模块进行误差剔除后粗检结果为影像I_R，在影像I_R中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为种子点；

第三子模块，用于若满足m>T且d<D则命令第四子模块工作，若满足m<T且k＝k1则命令第五子模块工作，否则将影像I_R中面积小于预设面积阈值的云区域清除，然后结束公子；

第四子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，命令第二子模块工作；

第五子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，命令第二子模块工作；

其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。

5.根据权利要求4所述多光谱遥感卫星影像自动云检测系统，其特征在于：云层粗提取模块中，所述基底阈值采用以下方式确定，

首先，计算Otsu阈值采用如下公式，

T＝MAX{U_t0×(V_to-V_tf)²+U_t1×(V_t1-V_tf)²}

其中，U_t0表示当阈值为t时，前景像元占总像元数比例，V_t0表示前景像元平均灰度，U_t1表示背景像元占总像元数比例，V_t1表示背景像元平均灰度，V_tf表示影像总平均灰度，V_tf＝U_t0×V_t0+U_t1×V_t1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t作为Otsu阈值T_J；

设云检测基底灰度阈值变化空间为[T_aT_b]，T_a、T_b分别为云检测基底灰度阈值变化下限和上限，按下式得到修正后的Otsu阈值作为基底阈值，

$T_{JF}=\left\{\begin{array}{cc}{T}_a,& if{T}_J<T_a\\ T_J,& if{T}_a\&le;T_J\&le;T_b\\ T_b,& if{T}_J>T_b\end{array}\right.$

其中T_J为Otsu阈值，T_JF为修正后的Otsu阈值。

6.根据权利要求4或5所述多光谱遥感卫星影像自动云检测系统，其特征在于：云层误差剔除模块中，

所述直方图均衡化的计算公式如下，

$S_k=\frac{\left({\mathrm{\&Sigma;}}_{j=0}^k{n}_j\right)}{N}\&times;L,k\&Element;\left(\begin{array}{cc}0& L\end{array}\right)$

公式中，S_k表示均衡化后的灰度值，表示灰度色阶0～k像素个数，k属于0到L；n_j表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶；

将直方图均衡化之后的强度灰度值进行双边滤波处理，包括将直方图均衡化之后强度S_k按下式得到滤波影像I’，

$I^{,}(i,j)=\frac{{\mathrm{\&Sigma;}}_{(i,j)\&Element;{\mathrm{\&Omega;}}_{x,y}}w(i,j)S_k(i,j)}{{\mathrm{\&Sigma;}}_{(i,j)\&Element;{\mathrm{\&Omega;}}_{x,y}}w(i,j)}$

w(i,j)＝W_s(i,j)×W_r(i,j)

$W_S(i,j)=e^{\frac{-\left|\right|(i,j)-(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_S}^2}}$

$W_r(i,j)=e^{\frac{-\left|\right|S_k(i,j)-S_k(i^{,},j^{,})|{|}^2}{2{{\mathrm{\&sigma;}}_r}^2}}$

D＝|S_k-I’|

其中，Ω_x,y表示中心点(x，y)的M×M大小的领域，对该领域内每一个像素点(i,j)均衡化后强度值S_k(i,j)，记w(i,j)为计算权重，W_s(i,j)为空间距离因子，W_r(i,j)为亮度相似度因子，σ_s和σ_r表示分别相应的高斯函数标准偏差，表示均衡化影像中强度最大值，e表示(i,j)为待滤波中心像素坐标，(i’,j’)为M*M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。