CN107066989A - 一种同步卫星遥感序列影像的积雪识别方法及系统 - Google Patents

Abstract

本发明公开了一种同步卫星遥感序列影像的积雪识别方法及系统。具体为，计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差，并基于光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像；计算动态时间规整值，并基于所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；基于积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像，充分利用了同步卫星遥感序列影像的多时相特性，无需短波红外波段，即可排除云的影响，实现同步卫星遥感影像中积雪的快速准确的识别。

Description

一种同步卫星遥感序列影像的积雪识别方法及系统

技术领域

本发明涉及遥感图像分类识别的领域，特别涉及一种同步卫星遥感序列影像的积雪识别方法及系统。

背景技术

雪的空间分布是山区与季节性积雪区水文、气象模型的重要输入因子，研究积雪分布及其特性，对于大气环流分析，气候演变，水资源利用以及灾害分析有着重要的意义。

李三妹等研究了国内静止气象卫星FY-2C资料的积雪判识原理，提出了基于阈值法的辅助因子函数积雪判识方法，较好地解决了FY-2C全圆盘范围内广大区域不同下垫面类型下的实时积雪监测问题。de Wildt等研究了基于Meteosat卫星数据的雪盖制图方法，利用静止卫星的高时相数据拓展了传统的积雪分类方法，更好地实现了云和雪的区分，得到了更加精准的积雪分布图。

地球同步轨道卫星具有高时间分辨率、观测范围广的优势，使得它能够实现对同一地区连续工作，能够更有效地实现实时地物类型-地物监测，一直是积雪监测十分重要的数据来源，根据这一特点地球同步轨道卫星在水资源利用、气候演变以及灾情监测等方面具有重要的应用价值。

2015年底，中国的GF-4(GaoFen-4)地球同步轨道光学卫星的成功发射，为区域雪盖制图提供了强大的数据支撑。GF-4卫星是我国首颗地球同步轨道高分辨率光学成像卫星，也是目前世界上空间分辨率最高、幅宽最大的地球同步轨道遥感卫星。由于GF-4卫星传感器不包含短波红外波段通道(1.55μm-1.65μm),NDSI方法也不再适用于GF-4卫星影像的积雪识别，因此积雪与云的区分成为新的难点。近年来使用缺少短波红外波段通道(1.55μm-1.65μm)的卫星影像，进行积雪识别的研究工作也有很多。

例如NOAA/AVHRR资料拥有5个观测通道包括了从可见光到远红外的光谱范围，但是这五个通道中并不包含有区分云雪的理想波段。殷青军等通过分析NOAA/AVHRR资料中云和雪的光谱特征，根据云和雪在热红外波段 AVHRR-3(3.55μm-3.93μm)和AVHRR-4(10.5μm-11.3μm)的亮温差，提出新的云、雪判别因子，对于区分青南地区的云、雪有较好的效果。但是在实际判识方面仍然有缺陷，如薄云的影响，较低云层的影响以及云层造成的阴影等。黄艳艳等针对国产卫星HJ-1B数据积雪像元的识别问题，根据积雪在HJ卫星CDD传感器波段上的反射率远高于其他地物，提出了一个仅用HJ-CDD数据作为数据源的积雪识别方法，该方法可以去除薄云的影响，但不能消除厚云的影响(黄艳艳,赵红莉,杨树红,2016.HJ-1B卫星遥感影像的积雪识别.测绘科学.41(8):129-133)。实际上，对于缺少短波红外通道的卫星而言，使用其影像作为唯一数据源提取积雪信息，很难消除厚云的影响。也有学者考虑使用多源数据进行积雪识别，曹云刚和刘闯指出，对于云层下面的信息，目前几乎所有的技术手段都难以将其完全恢复。比较好的途径是利用微波遥感数据提取积雪信息。但是这些方法有许多的弊端：一方面，预处理流程复杂，需要对不同的卫星数据进行几何校正，大气校正、裁剪等工作，并且两种数据的几何校正的精度很难保证极高的一致性。另一方面，由于空间分辨率，宽幅，扫描区域等均存在差异，获取同一区域、同一时间的两种卫星数据比较困难，所以该方法的实用性较低。

因此，如何根据同步卫星多时相遥感影像，排除云的干扰，实现卫星积雪的快速识别，成为亟待解决的技术问题。

发明内容

本发明的目的是，为了克服现有的GF-4卫星多时相遥感影像识别方法中云信息的干扰，提供了一种同步卫星遥感序列影像的积雪识别方法及系统。

为实现上述目的，本发明提供了如下方案：

一种同步卫星遥感序列影像的积雪识别方法，包括如下步骤：

根据同步卫星遥感影像库，选取无云信息影响的积雪影像，构建单时相积雪样本集和多时相积雪样本集；

在单时相积雪样本集中随机选取W个单时相积雪样本，计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差；

根据所述光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相 高亮地物影像；

将所述多时相高亮地物影像，按照时间序列进行叠加，生成多时相叠加影像；

根据所述多时相积雪样本集和所述多时相叠加影像，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值；

根据所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；

根据积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像。

可选的，所述构建单时相积雪样本集和多时相积雪样本集的步骤具体包括：

根据所述积雪样本影像中实际的地物类型定义A个积雪类别，所述积雪类别的编号为1～A，A为正整数，对于第i个积雪类别分别选取T_i个像元样本作为训练样本，构成单时相积雪样本，其中i＝1,2,…,A，T_i为正整数，且不大于第i类样本像元总数，由单时相积雪样本组成单时相积雪样本集；

基于同一地区，不同成像时间的M个单时相积雪样本，对于第i个积雪类别中的第j个单时相积雪样本，选取S_ij个像元，其中，j＝1,2，…，M，M为正整数，按成像时间将同一积雪类型，同一位置的像元的波段排列成向量，构成多时相积雪样本，其中S_ij为正整数，且不大于第i类样本像元总数的M倍，由多时相积雪样本组成多时相积雪样本集。

可选的，所述计算多时相待识别影像与单时相积雪样本之间的光谱角度的步骤具体包括：

按光谱角公式分别计算所述多时相待识别影像中的第i个像元与W个单时相样本中第i个单时相积雪样本之间的光谱角α_i,w，w＝1,2，…W，其中，X＝(a_i)_n为多时相待识别影像中的第i个像元光谱曲线矢量；Y＝(b_i,w)_n为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱曲线矢量；a_i为多时相待识别影像中第i个像元的光谱亮度， b_i,w为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度；n为多时相待识别影像所包含的波段数目；

将计算的W个光谱角，按升序排列，取第一个光谱角作为第i个像元与单时相积雪样本的光谱角度。

可选的，所述计算多时相待识别影像与单时相积雪样本之间的亮度差的步骤具体包括：

按亮度差值公式分别计算所述多时相待识别影像中的像元与W个单时相样本中第i个单时相积雪样本之间亮度差t_w，w＝1,2，…W，其中，a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度，n为多时相待识别影像所包含的波段数目；

将计算的W个亮度差，按升序排列，取第一个亮度差作为第i个像元与单时相积雪样本的亮度差。

可选的，所述根据光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像的步骤，具体包括：

利用一次最大类间方差算法，对光谱角度、亮度差值分别寻找到对应的初始阈值点和

分别构建光谱角度、亮度差的直方图，并得到直方图的中值点L_sam和L_lum，循环采用最大类间方差法得到光谱角度阈值I_SAM和亮度差阈值I_LUM；

根据阈值I_SAM和I_LUM，对于多时相待识别影像中光谱角度小于I_SAM并且亮度差小于I_LUM的像元，划归为疑似积雪类别，保留像元信息，否则划归为非积雪类别，将像元值置零，生成多时相高亮地物影像SM。

可选的，所述计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值的步骤具体包括：

利用动态时间规整算法计算多时相叠加影像的像元与多时相积雪样本之间的相似度；

对计算的相似度值进行升序排序，将第一个相似度值作为该像元与多时相积雪样本间的动态时间规整值。

可选的，所述对多时相高亮地物图像进行掩膜，生成积雪初步识别结果图像的步骤具体包括：

根据动态时间规整值，多次使用最大类间方差算法，自动寻找相似度的阈值ZDI；

判断多时相叠加影像的各像元的动态时间规整值是否大于相似度的阈值ZDI，若大于，则该像元为云，赋值为1，否则该像元为非云，赋值为0，得到云信息图像CM；

根据云信息图像CM，对多时相高亮地物图像进行掩膜，多时相高亮地物图像中对应位置的云信息图像CM的像元若为1，则将该位置的像元置零，否则，保留该位置的像元信息，生成积雪初步识别结果图像RM。

可选的，所述采用分类结果后处理的方法，得到积雪识别结果图像的步骤具体包括：

基于积雪初步识别结果图像，计算其每个像元的邻域范围内各像元类别出现的频数；选取频数最高的像元类别，标记为Lmax；

如果Lmax的频数为1，则不对当前像元像元进行处理，否则，赋予当前像元类别为Lmax；得到积雪识别结果图像RM。

发明还提供了一种同步卫星遥感序列影像的积雪识别系统，包括样本集构建模块，光谱角和亮度差计算模块，多时相待识别影像掩膜模块、影像叠加模块，动态时间规整模块，多时相高亮地物影像掩膜模块，分类结果后处理模块；

所述样本集构建模块，用于根据同步卫星遥感影像库，选取无云信息影响的积雪影像，构建单时相积雪样本集和多时相积雪样本集；

所述光谱角和亮度差计算模块，用于计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差；

所述多时相待识别影像掩膜模块，用于根据所述光谱角度和亮度差，对多 时相待识别影像进行掩膜，生成多时相高亮地物影像；

所述影像叠加模块，用于将多时相高亮地物影像，按照时间序列进行叠加，生成多时相叠加影像；

所述动态时间规整模块，用于根据多时相积雪样本集和多时相叠加影像，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值；

所述多时相高亮地物影像掩膜模块根据所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；

所述分类结果后处理模块，用于基于积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像。

可选的，所述样本集构建模块，包括单时相积雪样本集构建子模块和多时相积雪样本集构建子模块；

所述单时相积雪样本集构建子模块，用于根据所述积雪样本影像中实际的地物类型定义A个积雪类别，所述积雪类别的编号为1～A，A为正整数，对于第i个积雪类别分别选取T_i个像元样本作为训练样本，构成单时相积雪样本，其中i＝1,2,…,A，T_i为正整数，且不大于第i类样本像元总数，由单时相积雪样本组成单时相积雪样本集；

所述多时相积雪样本集构建子模块，用于基于同一地区，不同成像时间的M个单时相积雪样本，对于第i个积雪类别中的第j个单时相积雪样本，选取S_ij个像元，其中，j＝1,2，…，M，M为正整数，按成像时间将同一积雪类型，同一位置的像元的波段排列成向量，构成多时相积雪样本，其中S_ij为正整数，且不大于第i类样本像元总数的M倍，由多时相积雪样本组成多时相积雪样本集。

可选的，所述光谱角和亮度差计算模块，包括光谱角计算子模块和亮度差计算子模块；

所述光谱角计算子模块，用于按光谱角公式分别计算所述多时相待识别影像中的第i个像元与W个单时相样本中第i个单时 相积雪样本之间的光谱角α_i,w，w＝1,2，…W，其中，X＝(a_i)_n为多时相待识别影像中的第i个像元光谱曲线矢量；Y＝(b_i,w)_n为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱曲线矢量；a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度；n为多时相待识别影像所包含的波段数目；并将计算的W个光谱角，按升序排列，取第一个光谱角作为第i个像元与单时相积雪样本的光谱角度；

所述亮度差计算子模块，用于按亮度差值公式分别计算所述多时相待识别影像中的像元与W个单时相样本中第i个单时相积雪样本之间亮度差t_w，w＝1,2，…W，其中，a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度，n为多时相待识别影像所包含的波段数目；并将计算的W个亮度差，按升序排列，取第一个亮度差作为第i个像元与单时相积雪样本的亮度差。

可选的，所述多时相待识别影像掩膜模块包括光谱角度和亮度差阈值计算子模块和多时相待识别影像掩膜子模块；

所述光谱角度和亮度差阈值计算子模块，用于利用一次最大类间方差算法，对光谱角度、亮度差值分别寻找到对应的初始阈值点和并分别构建光谱角度、亮度差的直方图，并得到直方图的中值点L_sam和L_lum，循环采用最大类间方差法得到光谱角度阈值I_SAM和亮度差阈值I_LUM；

所述多时相待识别影像掩膜子模块，用于根据阈值I_SAM和I_LUM，对于多时相待识别影像中光谱角度小于I_SAM并且亮度差小于I_LUM的像元，划归为疑似积雪类别，保留像元信息，否则划归为非积雪类别，将像元值置零，生成多时相高亮地物影像SM。

可选的，所述动态时间规整模块包括相似度计算子模块和动态规整值计算子模块；

所述相似度计算子模块，用于利用动态时间规整算法计算多时相叠加影像的像元与多时相积雪样本之间的相似度；

所述动态规整值计算子模块，用于对计算的相似度值进行升序排序，将第一个相似度值作为该像元与多时相积雪样本间的动态时间规整值。

可选的，所述多时相高亮地物影像掩膜模块，包括相似度阈值计算子模块、云信息获取子模块和多时相高亮地物掩膜子模块；

所述相似度阈值计算子模块，用于根据动态时间规整值，多次使用最大类间方差算法，自动寻找相似度的阈值ZDI；

所述云信息获取子模块，用于判断多时相叠加影像的各像元的动态时间规整值是否大于相似度的阈值ZDI，若大于，则该像元为云，赋值为1，否则该像元为非云，赋值为0，得到云信息图像CM；

所述多时相高亮地物掩膜子模块，用于根据云信息图像CM，对多时相高亮地物图像进行掩膜，多时相高亮地物图像中对应位置的云信息图像CM的像元若为1，则将该位置的像元置零，否则，保留该位置的像元信息，生成积雪初步识别结果图像RM。

可选的，所述分类结果后处理模块，包括像元频数计算子模块和像元处理子模块，

所述像元频数计算子模块，用于基于积雪初步识别结果图像，计算其每个像元的邻域范围内各像元类别出现的频数；选取频数最高的像元类别，标记为Lmax；

像元处理子模块，用于根据Lmax像元的频数，对像元进行处理，如果Lmax的频数为1，则不对当前像元进行处理，否则，赋予当前像元类别为Lmax；得到积雪识别结果图像RM。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

发明公开了一种同步卫星遥感序列影像的积雪识别方法及系统，采用光谱角匹配法计算光谱角度，采用亮度差值方法计算亮度差；并基于光谱角度和亮 度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像；采用动态时间规整法得到动态时间规整值；并基于所述动态时间规整值，对多时相高亮地物图像进行掩膜，生成积雪初步识别结果图像；基于积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像，充分利用了同步卫星遥感序列影像的多时相特性，无需短波红外波段，即可排除云的影响，实现同步卫星遥感影像中积雪的快速准确的识别。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种同步卫星遥感序列影像的积雪识别方法的一个实施例的流程图。

图2为本发明提供的一种同步卫星遥感序列影像的积雪识别方法的另一个实施例的流程图。

图3为本发明提供的一种同步卫星遥感序列影像的积雪识别系统的结构框图。

具体实施方式

本发明的目的是提供一种面向群智多模态数据的处理方法及系统。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1所示，作为一种可实施方式，一种同步卫星遥感序列影像的积雪识别方法，包括如下步骤：

S1、根据同步卫星遥感影像库，选取无云信息影响的积雪影像，构建单时相积雪样本集和多时相积雪样本集；

S2、在单时相积雪样本集中随机选取W个单时相积雪样本，计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差；

S3、根据所述光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多 时相高亮地物影像；

S4、将所述多时相高亮地物影像，按照时间序列进行叠加，生成多时相叠加影像；

S5、根据所述多时相积雪样本集和所述多时相叠加影像，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值；

S6、根据所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；

S7、根据积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像。

如图2所示，可选的，作为一种可实施方式，步骤S2所述的计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差的方法为：根据单时相积雪样本集，基于时间和区域调用，采用光谱角匹配方法SAM和亮度差值方法计算多时相待识别图像中像元与单时相积雪样本之间的光谱角和亮度差，具体步骤包括：

按光谱角公式分别计算所述多时相待识别影像中的第i个像元与W个单时相样本中第i个单时相积雪样本之间的光谱角α_i,w，w＝1,2，…W，其中，X＝(a_i)_n为多时相待识别影像中的第i个像元光谱曲线矢量；Y＝(b_i,w)_n为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱曲线矢量；a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度；n为多时相待识别影像所包含的波段数目；

将计算的W个光谱角，按升序排列，取第一个光谱角作为第i个像元与单时相积雪样本的光谱角度。

按亮度差值公式分别计算所述多时相待识别影像中的像元与W个单时相样本中第i个单时相积雪样本之间亮度差t_w，w＝1,2，…W，其中，a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w第w个单时相积 雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度，n为多时相待识别影像所包含的波段数目；

将计算的W个亮度差，按升序排列，取第一个亮度差作为第i个像元与单时相积雪样本的亮度差。

步骤S3所述根据光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像的方法为多次利用改进的最大类间方差算法(OTSU算法)，自动寻找光谱阈值和亮度差阈值，并根据阈值掩膜多时相待识别图像得到多时相高亮地物图像，其中最大类间方差算法循环的停止条件为，寻找到阈值小于对应波段直方图统计的中值点，具体步骤包括：

利用一次最大类间方差算法，对光谱角度、亮度差值分别寻找到对应的初始阈值点和

分别构建光谱角度、亮度差的直方图，并得到直方图的中值点L_sam和L_lum，循环采用最大类间方差法得到光谱角度阈值I_SAM和亮度差阈值I_LUM；

根据阈值I_SAM和I_LUM，对于光谱角度小于I_SAM并且亮度差小于I_LUM的像元，划归为疑似积雪类别，保留像元信息，否则划归为非积雪类别，将像元值置零，生成多时相高亮地物影像SM。

所述最大类间方差算法为：设一副图像的灰度值为1～m级，在k处将其分成两组C₀＝{1～k}和C₁＝{k+1～m}，分别计算各组产生的概率w₀与w₁，各组的组内平均值μ₀与μ₁及整体图像的灰度平均值μ，两组间的方差用如下公式计算：

σ²(k)＝w₀(μ₀-μ)²+w₁(μ₁-μ)²＝w₀w₁(μ₁-μ₀)²＝[μ·w(k)-μ(k)]²/w(k)[1-w(k)]

从1～m之间改变k，求上式为最大值时的k，即求maxσ²(k)时的k*值便是所求的阈值；

所述得到光谱角度阈值和亮度差阈值的具体算法为：

若大于等于中值点L_sam(L_lum)，则转到步骤②计算若 小于中值点L_sam(L_lum)，则转到步骤③；

②令图像的灰度值范围为根据最大类间方差算法寻找新的阈值并转到步骤①；

③ISAM(ILUM)等于

步骤S5所述计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值的方法为利用动态规整方法(DTW)计算多时相叠加影像与多时相积雪样本之间的相似度，具体步骤包括：

利用动态时间规整算法计算多时相叠加影像的像元与多时相积雪样本之间的相似度；

对计算的相似度值进行升序排序，将第一个相似度值作为该像元与多时相积雪样本间的动态时间规整值。

所述动态时间规整算法为：假设有两个时间序列描述的特征A和B，A的长度为m，B的长度为n，其中：A＝{a1，a2，…，am}，B＝{b1，b2，…，bn}，于是DTW距离度量公式为：

DTW(A，B)＝min{cp(A，B)}

其中为了使用DTW让这两个时间序列非线性对齐，需要构造一个m×n的代价矩阵C＝[c(i，j)]，其中矩阵的第(i，j)元素为ai和aj的欧式距离c(i，j)。A与B之间规整路径p的累计代价函数可表示为

A与B之间的最优规整路径就是累计代价函数取得最小值时候的路径。

步骤S6所述对多时相高亮地物图像进行掩膜，生成积雪初步识别结果图像的方法为，根据多时相积雪样本集，基于时间和区域调用，采用最大类间方差法，对多时相高亮地物图像进行掩膜，得到积雪初步识别结果图像，具体步骤包括：

根据动态时间规整值，多次使用最大类间方差算法，自动寻找相似度的阈值ZDI；

判断多时相叠加影像的各像元的动态时间规整值是否大于相似度的阈值ZDI，若大于，则该像元为云，赋值为1，否则该像元为非云，赋值为0，得 到云信息图像CM；

根据云信息图像CM，对多时相高亮地物图像进行掩膜，多时相高亮地物图像中对应位置的云信息图像CM的像元若为1，则将该位置的像元置零，否则，保留该位置的像元信息，生成积雪初步识别结果图像RM。

所述相似度阈值寻找的具体步骤为：分别统计多时相动态时间规整值的直方图(直方图取值范围为：1～d)，并得到直方图的中值点ZDL；令i＝0,1,2….，于是寻找动态时间规整的阈值ZDI⁽ⁱ⁾：

①若ZDI⁽ⁱ⁾小于等于中值点ZDL，则转到步骤②计算ZDI⁽ⁱ⁺¹⁾；若ZDI⁽ⁱ⁾大于中值点ZDL，则转到步骤③；

②令图像的灰度值范围为ZD(i)～d，根据最大类间方差算法寻找新的阈值ZDI^{(i +1)}，并转到步骤①；

③ZDI等于ZDI⁽ⁱ⁾。

步骤S7所述采用分类结果后处理的方法，得到积雪识别结果图像的步骤具体包括：

基于积雪初步识别结果图像，计算其每个像元的邻域范围内各像元类别出现的频数；选取频数最高的像元类别，标记为Lmax；

如果Lmax的频数为1，则不对当前像元像元进行处理，否则，赋予当前像元类别为Lmax；得到积雪识别结果图像RM。

如图3所示，一种同步卫星遥感序列影像的积雪识别系统，包括样本集构建模块1，光谱角和亮度差计算模块2，多时相待识别影像掩膜模块3、影像叠加模块4，动态时间规整模块5，多时相高亮地物影像掩膜模块6，分类结果后处理模块7；

所述样本集构建模块1，用于根据同步卫星遥感影像库，选取无云信息影响的积雪影像，构建单时相积雪样本集和多时相积雪样本集；

所述光谱角和亮度差计算模块2，用于计算多时相待识别影像与单时相积 雪样本之间的光谱角度和亮度差；

所述多时相待识别影像掩膜模块3，用于根据所述光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像；

所述影像叠加模块4，用于将多时相高亮地物影像，按照时间序列进行叠加，生成多时相叠加影像；

所述动态时间规整模块5，用于根据多时相积雪样本集和多时相叠加影像，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值；

所述多时相高亮地物影像掩膜模块6，用于根据所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；

所述分类结果后处理模块7，用于基于积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像。

可选的，所述样本集构建模块1，包括单时相积雪样本集构建子模块和多时相积雪样本集构建子模块；

所述单时相积雪样本集构建模块，用于根据所述积雪样本影像中实际的地物类型定义A个积雪类别，所述积雪类别的编号为1～A，A为正整数，对于第i个积雪类别分别选取T_i个像元样本作为训练样本，构成单时相积雪样本，其中i＝1,2,…,A，T_i为正整数，且不大于第i类样本像元总数，由单时相积雪样本组成单时相积雪样本集；

所述多时相积雪样本集构建模块，用于基于同一地区，不同成像时间的M个单时相积雪样本，对于第i个积雪类别中的第j个单时相积雪样本，选取S_ij个像元，其中，j＝1,2，…，M，M为正整数，按成像时间将同一积雪类型，同一位置的像元的波段排列成向量，构成多时相积雪样本，其中S_ij为正整数，且不大于第i类样本像元总数的M倍，由多时相积雪样本组成多时相积雪样本集。

可选的，所述光谱角和亮度差计算模块2，包括光谱角计算子模块和亮度差计算子模块；

所述光谱角计算子模块，用于按光谱角公式分别计算所述多时相待识别影像中的第i个像元与W个单时相样本中第i个单时相积雪样本之间的光谱角α_i,w，w＝1,2，…W，其中，X＝(a_i)_n为多时相待识别影像中的第i个像元光谱曲线矢量；Y＝(b_i,w)_n为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱曲线矢量；a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度；n为多时相待识别影像所包含的波段数目；并将计算的W个光谱角，按升序排列，取第一个光谱角作为第i个像元与单时相积雪样本的光谱角度；

所述亮度差计算子模块，用于按亮度差值公式分别计算所述多时相待识别影像中的像元与W个单时相样本中第i个单时相积雪样本之间亮度差t_w，w＝1,2，…W，其中，a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度，n为多时相待识别影像所包含的波段数目；并将计算的W个亮度差，按升序排列，取第一个亮度差作为第i个像元与单时相积雪样本的亮度差。

可选的，所述多时相待识别影像掩膜模块3包括光谱角度和亮度差阈值计算子模块和多时相待识别影像掩膜子模块；

所述光谱角度和亮度差阈值计算子模块，用于利用一次最大类间方差算法，对光谱角度、亮度差值分别寻找到对应的初始阈值点和并分别构建光谱角度、亮度差的直方图，并得到直方图的中值点L_sam和L_lum，循环采用最大类间方差法得到光谱角度阈值I_SAM和亮度差阈值I_LUM；

所述多时相待识别影像掩膜子模块，用于根据阈值I_SAM和I_LUM，对于多时相待识别影像中光谱角度小于I_SAM并且亮度差小于I_LUM的像元，划归为疑似积雪类别，保留像元信息，否则划归为非积雪类别，将像元值置零，生成多 时相高亮地物影像SM。

可选的，所述动态时间规整模块5包括相似度计算子模块和动态规整值计算子模块；

所述相似度计算子模块，用于利用动态时间规整算法计算多时相叠加影像的像元与多时相积雪样本之间的相似度；

所述动态规整值计算子模块，用于对计算的相似度值进行升序排序，将第一个相似度值作为该像元与多时相积雪样本间的动态时间规整值。

可选的，所述多时相高亮地物影像掩膜模块6，包括相似度阈值计算子模块、云信息获取子模块和多时相高亮地物掩膜子模块；

所述相似度阈值计算子模块，用于根据动态时间规整值，多次使用最大类间方差算法，自动寻找相似度的阈值ZDI；

所述云信息获取子模块，用于判断多时相叠加影像的各像元的动态时间规整值是否大于相似度的阈值ZDI，若大于，则该像元为云，赋值为1，否则该像元为非云，赋值为0，得到云信息图像CM；

所述多时相高亮地物掩膜子模块，用于根据云信息图像CM，对多时相高亮地物图像进行掩膜，多时相高亮地物图像中对应位置的云信息图像CM的像元若为1，则将该位置的像元置零，否则，保留该位置的像元信息，生成积雪初步识别结果图像RM。

可选的，所述分类结果后处理模块7，包括像元频数计算子模块和像元处理子模块，

所述像元频数计算子模块，用于基于积雪初步识别结果图像，计算其每个像元的邻域范围内各像元类别出现的频数；选取频数最高的像元类别，标记为Lmax；

像元处理子模块，用于根据Lmax像元的频数，对像元进行处理，如果Lmax的频数为1，则不对当前像元进行处理，否则，赋予当前像元类别为Lmax；得到积雪识别结果图像RM。

本文中应用的实施例可以应用于GF-4卫星，排除云信息的影响，实现积雪的快速识别，但不限于GF-4卫星，也可以应用于其他的同步卫星，实现积雪的快速识别，本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

Claims (10)

1.一种同步卫星遥感序列影像的积雪识别方法，其特征在于，所述积雪识别方法包括

根据同步卫星遥感影像库，选取无云信息影响的积雪影像，构建单时相积雪样本集和多时相积雪样本集；

在单时相积雪样本集中随机选取W个单时相积雪样本，计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差；

根据所述光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像；

将所述多时相高亮地物影像，按照时间序列进行叠加，生成多时相叠加影像；

根据所述多时相积雪样本集和所述多时相叠加影像，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值；

根据所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；

根据积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像。

2.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，所述构建单时相积雪样本集和多时相积雪样本集的步骤具体包括：

根据所述积雪样本影像中实际的地物类型定义A个积雪类别，所述积雪类别的编号为1～A，A为正整数，对于第i个积雪类别分别选取T_i个像元样本作为训练样本，构成单时相积雪样本，其中i＝1,2,…,A，T_i为正整数，且不大于第i类样本像元总数，由单时相积雪样本组成单时相积雪样本集；

基于同一地区，不同成像时间的M个单时相积雪样本，对于第i个积雪类别中的第j个单时相积雪样本，选取S_ij个像元，其中，j＝1,2，…，M，M为正整数，按成像时间将同一积雪类型，同一位置的像元的波段排列成向量，构成多时相积雪样本，其中S_ij为正整数，且不大于第i类样本像元总数的M倍，由多时相积雪样本组成多时相积雪样本集。

3.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，所述计算多时相待识别影像与单时相积雪样本之间的光谱角度的步骤具体包括：

按光谱角公式分别计算所述多时相待识别影像中的第i个像元与W个单时相样本中第i个单时相积雪样本之间的光谱角α_i,w，w＝1,2，…W，其中，X＝(a_i)_n为多时相待识别影像中的第i个像元光谱曲线矢量；Y＝(b_i,w)_n为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱曲线矢量；a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度；n为多时相待识别影像所包含的波段数目；

将计算的W个光谱角，按升序排列，取第一个光谱角作为第i个像元与单时相积雪样本的光谱角度。

4.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，所述计算多时相待识别影像与单时相积雪样本之间的亮度差的步骤具体包括：

按亮度差值公式分别计算所述多时相待识别影像中的像元与W个单时相样本中第i个单时相积雪样本之间亮度差t_w，w＝1,2，…W，其中，a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度，n为多时相待识别影像所包含的波段数目；

将计算的W个亮度差，按升序排列，取第一个亮度差作为第i个像元与单时相积雪样本的亮度差。

5.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，根据所述光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像的步骤，具体包括：

利用一次最大类间方差算法，对光谱角度、亮度差值分别寻找到对应的初始阈值点和

分别构建光谱角度、亮度差的直方图，并得到直方图的中值点L_sam和L_lum，循环采用最大类间方差法得到光谱角度阈值I_SAM和亮度差阈值I_LUM；

根据阈值I_SAM和I_LUM，对于多时相待识别影像中光谱角度小于I_SAM并且亮度差小于I_LUM的像元，划归为疑似积雪类别，保留像元信息，否则划归为非积雪类别，将像元值置零，生成多时相高亮地物影像SM。

6.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值的步骤具体包括：

利用动态时间规整算法计算多时相叠加影像的像元与多时相积雪样本之间的相似度；

对计算的相似度值进行升序排序，将第一个相似度值作为该像元与多时相积雪样本间的动态时间规整值。

7.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，所述对多时相高亮地物图像进行掩膜，生成积雪初步识别结果图像的步骤具体包括：

根据动态时间规整值，多次使用最大类间方差算法，自动寻找相似度的阈值ZDI；

判断多时相叠加影像的各像元的动态时间规整值是否大于相似度的阈值ZDI，若大于，则该像元为云，赋值为1，否则该像元为非云，赋值为0，得到云信息图像CM；

根据云信息图像CM，对多时相高亮地物图像进行掩膜，多时相高亮地物图像中对应位置的云信息图像CM的像元若为1，则将该位置的像元置零，否则，保留该位置的像元信息，生成积雪初步识别结果图像RM。

8.根据权利要求1所述的一种同步卫星遥感序列影像的积雪识别方法，其特征在于，所述采用分类结果后处理的方法，得到积雪识别结果图像的步骤具体包括：

基于积雪初步识别结果图像，计算其每个像元的邻域范围内各像元类别出现的频数；选取频数最高的像元类别，标记为Lmax；

如果Lmax的频数为1，则不对当前像元像元进行处理，否则，赋予当前像元类别为Lmax；得到积雪识别结果图像RM。

9.一种同步卫星遥感序列影像的积雪识别系统，其特征在于，包括样本集构建模块，光谱角和亮度差计算模块，多时相待识别影像掩膜模块、影像叠加模块，动态时间规整模块，多时相高亮地物影像掩膜模块，分类结果后处理模块；

所述样本集构建模块，用于根据同步卫星遥感影像库，选取无云信息影响的积雪影像，构建单时相积雪样本集和多时相积雪样本集；

所述光谱角和亮度差计算模块，用于计算多时相待识别影像与单时相积雪样本之间的光谱角度和亮度差；

所述多时相待识别影像掩膜模块，用于根据所述光谱角度和亮度差，对多时相待识别影像进行掩膜，生成多时相高亮地物影像；

所述影像叠加模块，用于将多时相高亮地物影像，按照时间序列进行叠加，生成多时相叠加影像；

所述动态时间规整模块，用于根据多时相积雪样本集和多时相叠加影像，计算多时相叠加影像与多时相积雪样本之间的相似度，得到动态时间规整值；

所述多时相高亮地物影像掩膜模块根据所述动态时间规整值，对多时相高亮地物影像进行掩膜，生成积雪初步识别结果图像；

所述分类结果后处理模块，用于基于积雪初步识别结果图像，采用分类结果后处理的方法，得到积雪识别结果图像。

10.根据权利要求9所述的一种同步卫星遥感序列影像的积雪识别系统，其特征在于，所述光谱角和亮度差计算模块，包括光谱角计算子模块和亮度差计算子模块；

所述光谱角计算子模块，用于按光谱角公式分别计算所述多时相待识别影像中的第i个像元与W个单时相样本中第i个单时相积雪样本之间的光谱角α_i,w，w＝1,2，…W，其中，X＝(a_i)_n为多时相待识别影像中的第i个像元光谱曲线矢量；Y＝(b_i,w)_n为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱曲线矢量；a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w为第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度；n为多时相待识别影像所包含的波段数目；并将计算的W个光谱角，按升序排列，取第一个光谱角作为第i个像元与单时相积雪样本的光谱角度；

所述亮度差计算子模块，按亮度差值公式分别计算所述多时相待识别影像中的像元与W个单时相样本中第i个单时相积雪样本之间亮度差t_w，w＝1,2，…W，其中，a_i为多时相待识别影像中第i个像元的光谱亮度，b_i,w第w个单时相积雪样本中与多时相待识别影像中的第i个像元对应的像元的光谱亮度，n为多时相待识别影像所包含的波段数目；并将计算的W个亮度差，按升序排列，取第一个亮度差作为第i个像元与单时相积雪样本的亮度差。