[0048] 第四子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，命令第二子模块工作；
[0049] 第五子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，命令第二子模块工作；
[0050] 其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。
[0051] 而且，云层粗提取模块中，所述基底阈值采用以下方式确定，
[0052] 首先，计算Otsu阈值采用如下公式，
[0053] T＝MAX{Ut0×(Vto-Vtf)2+Ut1×(Vt1-Vtf)2}
[0054] 其中，Ut0表示当阈值为t时，前景像元占总像元数比例，Vt0表示前景像元平均灰度，Ut1表示背景像元占总像元数比例，Vt1表示背景像元平均灰度，Vtf表示影像总平均灰度，Vtf＝Ut0×Vt0+Ut1×Vt1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t作为Otsu阈值TJ；
[0055] 设云检测基底灰度阈值变化空间为[Ta Tb]，Ta、Tb分别为云检测基底灰度阈值变化下限和上限，按下式得到修正后的Otsu阈值作为基底阈值，
[0056]
[0057] 其中TJ为Otsu阈值，TJF为修正后的Otsu阈值。
[0058] 而且，云层误差剔除模块中，
[0059] 所述直方图均衡化的计算公式如下，
[0060]
[0061] 公式中，Sk表示均衡化后的灰度值， 表示灰度色阶0～k像素个数，k属于0到L；nj表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶；
[0062] 将直方图均衡化之后的强度灰度值进行双边滤波处理，包括将直方图均衡化之后强度Sk按下式得到滤波影像I’，
[0063]
[0064] w(i,j)＝Ws(i,j)×Wr(i,j)
[0065]
[0066]
[0067] D＝|Sk-I’|
[0068] 其中，Ωx,y表示中心点(x，y)的M×M大小的领域，对该领域内每一个像素点(i,j)均衡化后强度值Sk(i,j)，记w(i,j)为计算权重，Ws(i,j)为空间距离因子，Wr(i,j)为亮度相似度因子，σs和σr表示分别相应的高斯函数标准偏差， 表示均衡化影像中强度最大值，e表示(i,j)为待滤波中心像素坐标，(i’,j’)为M*M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。
[0069] 本发明所采用多光谱阈值与纹理信息相结合的方法进行云检测，首先采用改进的颜色模型转换方法，将影像由RGB转换至HIS颜色空间。利用影像强度信息与饱和度信息结合Otsu阈值分割生成基底图，并结合影像近红外信息与色调信息对基底图进行优化，生成修正图。然后利用直方图均衡化与双边滤波结合二维Otsu分割方法提取地物的纹理信息，对修正图进行粗差剔除，生成云种子图。最后将强度信息作为向导结合种子图信息进行云精确提取。本发明主要应用于卫星影像云检测，保证卫星影像产品生产等后续工作的顺利进行、大范围高质量无云遥感影像合成等。本发明技术方案具有以下优点：
[0070] 1)利用改进的HIS模型，将原影像变换到HIS色彩空间，充分利用了云的光谱特性，大大降低了薄云的检测难度。
[0071] 2)利用直方图均衡化和双边滤波相结合的方法进行纹理信息提取，充分体现了云层粗纹理特征，能够将云与高相似度下垫面有效区分。
[0072] 3)二维Otsu阈值分割法，提高了阈值的适用性，并提高了纹理检测精度。
[0073] 4)边缘种子膨胀，大大提高了云层检测精度。
附图说明
[0074] 图1是本发明实施例的流程图。
[0075] 图2是本发明实施例的云层基底像素值统计示意图，其中图2(a)为灰度直方图，图2(b)为经验累积分布图。
[0076] 图3是本发明实施例的二维Otsu阈值生成原理示意图，其中图3(a)为一维Otsu阈值示意图，图3(b)为二维Otsu阈值示意图。
[0077] 图4是本发明实施例的膨胀操作原理示意图。
[0078] 图5是本发明实施例的云层精提取流程图。
具体实施方式
[0079] 以下结合实施例和附图详细说明本发明技术方案。
[0080] 本发明提供一种快速、准确、全自动化的基于阈值分割与纹理信息的多光谱遥感卫星影像自动云检测方法，主要利用影像的光谱以及纹理信息进行云提取，首先利用光谱信息进行云结果初步检测，然后利用纹理信息对初检结果中的误差进行剔除，最后回头利用光谱信息对云层进行精确检测。输入数据为能提供红、绿、蓝、近红外波段数据的各类高分辨遥感卫星影像。参见图1，实施例所提供具体实现方法包含以下步骤：
[0081] 步骤1：数据准备。获取待进行云检测的卫星影像，该影像中必须包含红、绿、蓝、近红外四个波段的数据。
[0082] 本方法适用于所有能够提供红、绿、蓝、近红外四个波段的高分辨率遥感卫星影像数据。研究发现，RGB彩色模型对厚云检测效果较好，但难以识别薄云，并且人的视觉对亮度的敏感程度远强于对颜色浓淡的敏感程度，因此HIS色彩空间比RGB色彩空间更符合人的视觉特性。因此，后续在步骤2通过改进的HIS模型仿真人眼感知色彩的特性，充分利用云的光谱特性，大大降低了薄云的检测难度。
[0083] 步骤2：云层粗提取。将原影像变换到HIS色彩空间，将强度信息I、饱和度信息S线性拉伸至[0 1]。然后计算得到基底影像。具体实施时，近红外阈值、色调阈值可根据统计结果预设，例如从高分一号、资源三号等卫星影像中，选取500个云对象(包含薄云、厚云、卷云等)，统计发现，云近红外值均大于350，色调值均小于120。通过带限定条件的Otsu阈值(作为基底阈值)结合近红外、色调信息的限定条件对基底影像进行分割，即基底像素值大于基底阈值(J大于TJF)并且近红外像素值大于近红外阈值并且色调像素值小于色调阈值的像元即为粗检云对象，而不符合任意某一条件的像元，即为非云对象，以此得到云层粗检结果，即对基底图去除无云区域的修正图。
[0084] 本发明涉及的RGB转HIS色彩通道计算公式，如以下所示：
[0085]
[0086] I＝1/3(R’+G’+B’)
[0087]
[0088]
[0089]
[0090] 式中R、G、B为原始影像色彩灰度值，R’为R、G、B中最小值，G’为R、G、B中间值，B’为R、G、B中最大值，min(R’,G’,B’)为R’,G’,B’中最小值，H0为影像色调值的初始值，H为影像色调值，S为影像饱和度信息，I为影像强度信息(即亮度)。
[0091] 本发明使用影像强度、饱和度、色调、近红外信息进行云层粗检测，其策略如下：
[0092] 将原影像变换到HIS色彩空间，并将I、S线性拉伸至[0 1](即变换到0-1之间)，并利用下述公式计算得到基底影像。
[0093]
[0094] 在基于像素进行计算时，公式中影像强度信息I’代表原影像中对应像素的强度值线性拉伸后结果，影像饱和度信息S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，具体实施时可以由本领域技术人员预设为大于0的值，实施例中取值为1.0，J代表生成的基底影像对应像素值，大小属于0-255。
[0095] 利用带限定条件的Otsu阈值法，对基底影像分割，生成云层粗检测结果。计算公式如下：
[0096] T＝MAX{Ut0×(Vto-Vtf)2+Ut1×(Vt1-Vtf)2}
[0097] 其中Ut0表示当阈值为t时，前景像元占总像元数比例，Vt0表示前景像元平均灰度，Ut1表示背景像元占总像元数比例，Vt1表示背景像元平均灰度，Vtf表示影像总平均灰度，Vtf＝Ut0×Vt0+Ut1×Vt1。从最小灰度值到最大灰度值遍历t，当t使得类间方差T最大时，此时t即为Otsu阈值TJ。Otsu阈值大体可将云层与下垫面分割开，但往往有分割阈值过大，薄云漏检；分割阈值过小误检率过高等现象存在，通过带约束的Otsu阈值能有效解决这一问题。
[0098] 根据云检测基底灰度阈值变化空间，可以对阈值进行修正，设云检测基底灰度阈值变化空间为[Ta Tb]，Ta、Tb分别为云检测基底灰度阈值变化下限和上限，按下式得到修正阈值，
[0099]
[0100] 具体实施时，云检测基底灰度阈值变化下限和上限可根据统计结果预设，例如从高分一号、资源三号等卫星影像中，选取500个云对象(包含薄云、厚云、卷云等)，统计其灰度直方图和经验累积分布图如图2(a)、图2(b)所示。
[0101] 根据图2分析，云层基底像素值大体分布在[80255]之间，累积频率在130开始呈现快速上升的趋势，表征当某基底像素灰度值为130时，为云的概率极大，此时如果计算Otsu阈值大于130，那么漏掉薄云的可能性将随之迅速增大。因此，实验中设定云检测基底灰度阈值变化空间为[80130]，即按下式得到修正阈值。
[0102]
[0103] 其中TJ为Otsu阈值，TJF为修正后的Otsu阈值。
[0104] 步骤3：云层误差剔除。在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方法进行纹理信息提取。然后利用二维Otsu阈值Td(即在一维Otsu阈值基础上，将满足阈值条件的像素作为下一轮直方图统计数据来源，再次计算Otsu阈值)对纹理信息图进行分割，利用分割后的细节图对步骤2所得粗检结果进行误差剔除，即根据细节图进行粗差剔除得到云种子图。
[0105] 本发明首先使用直方图均衡化突出隐含有纹理细节的图像。通过控制图像灰度级概率密度函数，改变图像的灰度层次。均衡化之后，山区灰度跳变幅度变大，城区呈现明显颗粒感，突出了建筑物等细节。而云层只是对比度增强，纹理却没有变细，灰度空间分布仍保持较好的连续性。
[0106] 直方图均衡化计算公式如下：
[0107]
[0108] 公式中，Sk表示均衡化后的灰度值， 表示灰度色阶0～k像素个数，k属于0到L；nj表示原图灰度色阶为j的像素个数，j的取值范围为0～k，N是图像像素总数，L表示均衡化之后的图像色阶。
[0109] 然后将直方图均衡化之后的强度灰度值进行双边滤波处理，双边滤波对像素同时使用空间和边界滤波，在对图像进行光滑的同时能保持图像边缘特征。对于细节信息较少的卫星影像，通常需要多次选权双边滤波提取地物纹理，但对计算机内存及计算速度是一个极大考验。文中直方图均衡化处理，有效强化了地物纹理，仅需一次双边滤波就能有效提取地物纹理。将直方图均衡化之后强度Sk按下式得到滤波影像I’：
[0110]
[0111] w(i,j)＝Ws(i,j)×Wr(i,j)
[0112]
[0113]
[0114] D＝|Sk-I’|
[0115] 式中Ωx,y表示中心点(x，y)的M×M(M为奇数)大小的领域。对该领域内每一个像素点(i,j)均衡化后强度值Sk(i,j)，记w(i,j)为计算权重，Ws(i,j)为空间距离因子，Wr(i,j)为亮度相似度因子，σs和σr表示分别相应的高斯函数标准偏差，具体实施时，本领域技术人员可自行预设取值，实施例中σs取值为2，σr取值为 表示均衡化影像中强度最大值，e表示数学常量2.718281，是自然对数函数的底数，(i,j)为待滤波中心像素坐标，(i’,j’)为M×M窗口领域其他像素坐标，I’(i，j)为双边滤波后的中心像素强度值，D为细节图像素灰度值。
[0116] 采用二维Otsu阈值对细节图进行二值分割，即在一维Otsu阈值基础上，将满足阈值条件的像素作为下一轮直方图统计数据来源，代入Otsu阈值计算公式，再次计算Otsu阈值，过程如图3所示，如图3(a)中一维Otsu阈值Otsu_1＝22，如图3(b)中二维Otsu阈值Otsu_2＝7。具体实施时，可按照获取TJ方式执行一遍，将高于阈值的像素剔除，低于阈值条件的像素保留。根据保留下来的像素及相应统计直方图数据代入Otsu阈值计算公式，再次计算TJ。例如Ut0表示当阈值为t时前景像元占总像元数比例，再次计算时采用保留的像素相应统计直方图中所有大于t灰度的直方图部分的比例。
[0117] 二维Otsu阈值分割在一定程度上提高阈值的适用性，并提高了检测精度。随后，为剔除弱噪声，对分割结果进行一次3*3窗口的膨胀操作获取二值细节图。用二值细节图对粗检结果进行粗差剔除，计算公式如下：
[0118] IR＝IX∩ID
[0119] 式中IR表示粗差剔除后粗检结果影像，即为云种子图，IX表示粗检结果影像，ID表示二值细节图。
[0120] 其中，3*3窗口的膨胀操作原理如图4所示，将大小为3*3的模板与原始图卷积，若判断像素与其八领域均为检测目标则保留，否则剔除。
[0121] 步骤4：云层精提取。将步骤2所得原始影像强度信息作为向导图，在步骤3进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，得到最终的云检结果。当k＝0.8‰时，能够将漏检厚云提取出来。而先将厚云向周围薄云过渡，调整灰度阈值后，能够有效将厚云周边薄云提取出来。
[0122] 将步骤2所得RGB转HIS所得影像强度信息I作为向导图，在粗差剔除后粗检结果影像IR上利用边缘种子膨胀的方式进行精确提取，实施例的云层精确提取流程如图5所示，步骤如下：
[0123] 步骤4.1：初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为0.8‰，设定新增面积阈值T＝200，迭代次数阈值D＝10；具体实施时，k、T、D可由用户自行预设，例如实施例采用的实验统计最佳设定值；
[0124] 步骤4.2：粗差剔除后粗检结果影像IR中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为待处理的种子点；
[0125] 步骤4.3：若满足条件一：m>T且d
[0126] 步骤4.4：在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m、d(d＝d+1)值，返回步骤4.2；
[0127] 步骤4.5：在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k为10％作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值。所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k为1.2‰，返回步骤4.2。
[0128] 在步骤4.5中，通过调整k值，首先k取极小值k1＝0.8‰，将厚云周边的漏检的同类厚云搜索出来；然后将k调整为极大值k2＝10％，等同于强行将检测结果向外扩充一圈，将厚云向周边的漏检薄云过度；最后将k调整为较小值k3＝1.2‰，再次搜索，可将薄云周边漏检的同类薄云搜索出来。具体实施时，k1、k2、k3的取值可由本领域技术人员预设为经验值，建议k1,k3一般不超过1.5％,k2不超过15％。
[0129] 从理论上分析，该方法以光谱结合纹理的方法进行云检测。但不同于传统阈值法的是：其一：传统光谱结合阈值法，多采用RGB彩色模型，而本方法通过改进的HIS模型，能够有效降低薄云的检测难度；其二：该方法不拘泥于阈值的选取，通过带限定条件的Otsu阈值可保证厚、薄云被提取的同时，减少同类地物的干扰；其三：阈值分割过程中，引入近红外波段信息，在区分云与高亮的河流、积雪等地物时取得理想的效果。亦不同于传统纹理分析法的是：采用直方图均衡化与双边滤波相结合提取影像纹理信息的方法，有效的弥补了传统分形维数和角二阶矩计算精度低、效率低下的不足，且对厚云、薄云均能取得理想效果。
[0130] 利用本发明进行试验，包括全自动的对资源三号、高分一号、高分二号等高分辨率卫星影像进行高效率云检测，可验证分析该方法有以下优点：
[0131] 无需任何先验知识，全自动化处理；
[0132] 算法复杂度适中，计算量小，计算速度快；
[0133] 检测精度高，可在高精确度检测出厚云、薄云的前提下，有效排除高亮度水域、房屋、裸地等干扰；
[0134] 算法适用性广，可用于能提供红、绿、蓝、近红外波段数据的各类卫星影像的云检测。
[0135] 本发明技术方案可采用计算机软件方式支持自动运行流程，具体实施时也可采用模块化方式提供相应系统。本发明实施例提供一种多光谱遥感卫星影像自动云检测系统，包括以下模块：
[0136] 数据准备模块，用于获取待进行云检测的卫星影像，卫星影像中包含红、绿、蓝、近红外四个波段的数据；
[0137] 云层粗提取模块，用于将原影像变换到HIS色彩空间，将强度信息和饱和度信息分别线性拉伸至[0 1]，然后计算得到基底影像，判断基底像素值大于基底阈值且近红外像素值大于近红外阈值且色调像素值小于色调阈值的像元为粗检云对象，否则为非云对象，得到云层粗检结果；
[0138] 将强度信息和饱和度信息线性拉伸至[0 1]后，利用下述公式计算得到基底影像，[0139]
[0140] 其中，I’代表原影像中对应像素的强度值线性拉伸后结果，S’代表原影像对应像素的饱和度值线性拉伸后结果，τ为缓冲系数，J代表生成的基底影像对应像素值；
[0141] 云层误差剔除模块，用于在影像强度信息通道上，使用直方图均衡化和双边滤波相结合的方式进行纹理信息提取，利用二维Otsu阈值对纹理信息图进行分割，利用分割后所得二值细节图对粗检结果进行误差剔除，
[0142] 云层精提取模块，用于将云层粗提取模块变换所得原影像的强度信息作为向导图，在云层误差剔除模块进行误差剔除后粗检结果的基础上利用边缘种子膨胀的方式对云层进行精确提取，包括以下子模块，
[0143] 第一子模块，用于初始化新增云像元数m为正无穷，迭代次数d＝0，设定灰度阈值k为k1，设定新增面积阈值，迭代次数阈值D；
[0144] 第二子模块，用于设在云层误差剔除模块进行误差剔除后粗检结果为影像IR，在影像IR中检索每一云像元，判断其八邻域是否均为云像元，若是则不处理，否则判定为种子点；
[0145] 第三子模块，用于若满足m>T且d
[0146] 第四子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，若种子点与某一邻域点的灰度差值小于阈值k×I(i,j)，则将该邻域点判为新增云对象，I(i,j)为种子点原始影像强度值；所有种子点处理完成后，统计并更新m值，更新d＝d+1，命令第二子模块工作；
[0147] 第五子模块，用于在向导图上，分别计算各种子点与其八个邻域点的灰度差值，设定k＝k2作为判断依据，若种子点与某一邻域灰度差值小于k×I(i,j)，则将该邻域点判为新增薄云像元，I(i,j)为种子点原始影像强度值，所有种子点处理完成后，更新m值为正无穷，调整灰度阈值k＝k3，命令第二子模块工作；
[0148] 其中，k1为用于搜索厚云周边漏检的同类厚云的预设极小值，k2为用于搜索厚云周边的漏检薄云的预设极大值，k3为搜索薄云周边漏检的同类薄云的预设较小值。
[0149] 各模块具体实现可参加相应步骤，本发明不予赘述。
[0150] 本文中所描述的具体实施例仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例作各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或超越所附权利要求书所定义的范围。
