CN106127717B - 一种基于modis逐日积雪覆盖图像的去云方法及装置 - Google Patents

Abstract

本发明公开了一种基于MODIS的逐日积雪覆盖图像的去云方法及装置，方法包括：合成MODIS的上下午星逐日积雪覆盖图像；将连续三天的积雪覆盖图像进行合成，分别对合成图像进行基于“永久性”陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，将连续八天的合成图像进行“最大陆地”和“最大积雪”合成处理并根据处理结果进行云覆盖像素的积雪判别后，根据拟合预期雪线方法对进行去云处理，得到无云积雪图像。本发明通过将MODIS上下午卫星逐日积雪图像数据进行合成，并根据一系列的阈值等多步骤判别和拟合方法对合成后的积雪图像进行去云操作，能够去除逐日积雪图像中影响动态监测的云，满足监测者的要求。

Description

一种基于MODIS逐日积雪覆盖图像的去云方法及装置

技术领域

本发明涉及卫星积雪遥感产品处理技术领域，具体涉及一种基于MODIS逐日积雪覆盖图像的去云方法及装置。

背景技术

积雪是冰冻圈的重要组成部分，是地球表面最为活跃的自然要素之一，对区域甚至全球的气候变化、能量平衡以及水循环都有着深刻影响。青藏高原素有“世界屋脊”之称，是中纬度地区海拔最高的区域，也是“高亚洲地区”的核心区域。高原积雪对于高亚洲地区的区域气候和能量循环起着重要的反馈和调节作用，其消融过程也影响着融雪性河流流量的变化，对于下游地区的水、生态和灾害具有重要的作用，对社会和经济产生重要的影响。

青藏高原地形的复杂引起热量水汽不均，因此积雪分布也是极其不均匀的。积雪主要集中在西部的喀喇昆仑山山脉，南部的喜马拉雅山，以及高原北部的天山，东南部的念青唐古拉山脉附近。高原中部腹地包括柴达木盆地，以及北部的临近高原的塔里木盆地，由于受四周高山地形的影响，水汽输送较少，是积雪的小值区域。青藏高原地区独特复杂的地形以及积雪破碎的情况，使得积雪遥感动态监测面临着不小的挑战。

在青藏高原积雪研究中，由于地面观测站点较少，且代表性和观测时间长度有限，可通过卫星存档数据开展有效地积雪反演和制图研究，目前可以使用的雪盖估算遥感数据，可以追溯近40年前状态，是对地面台站的最有效的补充，并逐渐成为高原区常规性积雪监测的主要手段。利用光学遥感进行积雪面积监测已有50多年的历史，常见的积雪产品有Landsat、SPOT、AVHRR、和MODIS等。其中由Terra和Aqua卫星携带的中分辨率成像光谱仪MODIS，以其高时空分辨率、长运行时间及较好的稳定性得到了最为广泛的应用。青藏高原每日MODIS积雪数据(MOD10A1、MYD10A1)可以在日时间尺度上实现积雪覆盖面积的估算，为模型提供高时间分辨率的下垫面初始场，在一定程度上还可以帮助大尺度被动微波遥感提高积雪的判识能力。然而，由于积雪和云的反射光谱特征相似，MODIS日积雪产品受云污染严重，使积雪制图难以达到较高的精度。针对这一问题，许多学者开展了大量去除云污染的研究，提出多种去云算法，以期实现每日无云积雪产品的开发。

目前，卫星遥感技术由于重访其周期短(半月至日访问周期)、覆盖范围大等优点逐渐成为积雪参数监测的主要方式。当前，常用的积雪遥感数据有Landsat和SPOT、AVHRR、VEGETATION、MODIS等，可用于光学积雪覆盖度和面积产品，SMMR、SSM/I和AMSR-E等微波遥感数据，用于雪深监测。其中MODIS逐日积雪监测数据(MOD10A1和MYD10A1)以较高的空间分辨率和时间分辨率为特征，在北半球大面积积雪制图方面，其精度可达94％，可以满足全球气候和环境变化研究，且卫星运行较为稳定，应用最为广泛。然而由于积雪和云的反射特性相似，使得光学遥感产品受天气影响严重，在日时间尺度上，云污染限制了积雪产品的应用。

现有的卫星数据去云方法的去云效果较差，无法满足逐日时间尺度的动态监测要求。

发明内容

由于现有的卫星数据去云方法的去云效果较差，无法满足逐日时间尺度的动态监测的要求的问题，本发明提出一种基于MODIS的逐日积雪覆盖图像的去云方法及装置。

第一方面，本发明提出一种基于MODIS的逐日积雪覆盖图像的去云方法，包括：

获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；

根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；

将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于“永久性”陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；

根据第二合成优先级，将连续八天的所述第三积雪覆盖合成图像进行“最大陆地”和“最大积雪”合成处理，得到八天合成处理结果，并根据所述八天处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；

对所述第四积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

优选地，还包括：

根据预设积雪覆盖精度指标，对所述MODIS无云积雪图像进行验证；

其中，所述预设积雪覆盖精度指标包括准确度、召回率、精度和f值。

优选地，所述基于“永久性”陆地和积雪覆盖的判别处理进一步包括：

将所述第二积雪覆盖合成图像中，高程大于第一预设值的区域对应的像素点调整为积雪像素；

若当前像素为云像素，所述云像素的高程大于第二预设值小于等于所述第一预设值，且所述云像素的云天数与积雪天数之和大于总天数与第三预设值的乘积，则将所述云像素调整为积雪像素；

若所述云像素的云天数与陆地天数之和与所述总天数相同，且所述云天数小于所述总天数与第四预设值的乘积，则将所述云像素调整为陆地像素。

优选地，所述临近四像元处理进一步包括：

若所述云像素相邻的四个像素中至少三个为同类型非云像素，则将当前云像素调整为所述同类型非云像素；

其中，像素的数据类型包括：云、积雪、陆地、湖泊和湖冰；

所述同类型非云像素包括积雪像素，陆地像素、湖泊像素和湖冰像素。

优选地，所述高程滤波处理进一步包括：

若当前像素为积雪像素，所述积雪像素相邻的八个像素中有云像素，且所述云像素对应的高程大于所述积雪像素的高程，则将所述云像素调整为积雪像素。

优选地，所述将高山冰川覆盖区的湖水误判为阴影区处理进一步包括：

根据现有湖泊辅助数据像素的矢量边界，将所述矢量边界之外的湖泊像素和湖冰像素均调整为云像素。

优选地，所述拟合预期雪线方法进一步包括：

根据云像素的坡度、坡向、经度和纬度，得到云像素的高程曲线，并根据云的实际高程，将小于实际高程的云像素调整为积雪像素。

第二方面，本发明还提出一种基于MODIS的逐日积雪覆盖图像的去云装置，包括：

上下午合成模块，用于获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；

三天合成模块，用于根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；

处理后合成模块，用于将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于“永久性”陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；

八天合成模块，用于根据第二合成优先级，将连续八天的所述第三积雪覆盖合成图像进行“最大陆地”和“最大积雪”合成处理，得到八天合成处理结果，并根据所述八天处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；

雪线拟合模块，用于对所述第四积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

优选地，还包括：

图像验证模块，用于根据预设积雪覆盖精度指标，对所述MODIS无云积雪图像进行验证；

其中，所述预设积雪覆盖精度指标包括准确度、召回率、精度和f值。

优选地，所述处理后合成模块中积雪覆盖的判别处理进一步包括：

第一调整单元，用于将所述第二积雪覆盖合成图像中，高程大于第一预设值的区域对应的像素点调整为积雪像素；

第二调整单元，用于若当前像素为云像素，所述云像素的高程大于第二预设值小于等于所述第一预设值，且所述云像素的云天数与积雪天数之和大于总天数与第三预设值的乘积，则将所述云像素调整为积雪像素；

第三调整单元，用于若所述云像素的云天数与陆地天数之和与所述总天数相同，且所述云天数小于所述总天数与第四预设值的乘积，则将所述云像素调整为陆地像素。

由上述技术方案可知，本发明通过将多个卫星数据进行合成，并根据阈值对合成后的卫星图像进行多步去云操作，能够完全去除图像中影响动态监测的云，满足监测者的要求。

附图说明

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

图1为本发明一实施例提供的一种基于MODIS的逐日积雪覆盖图像的去云方法的流程示意图；

图2为本发明一实施例提供的MODIS逐日无云积雪产品合成算法流程图；

图3为本发明一实施例提供的MODIS逐日无云积雪产品永久积雪、陆地判别算法流程图；

图4为本发明一实施例提供的2010年10月1日-2011年4月30日“永久积雪陆地”分布；

图5为本发明一实施例提供的MODIS逐日无云积雪产品临近四像元算法流程图；

图6为本发明一实施例提供的2011年1月9日-2011年1月17日最小积雪面积、最小陆地面积及积雪变化范围；

图7为本发明一实施例提供的MODIS逐日无云积雪产品拟合雪线算法流程图；

图8为本发明一实施例提供的青藏高原分区图；

图9为本发明一实施例提供的2009-2011年积雪季过程产品云量分析图；

图10为本发明一实施例提供的一种基于MODIS的逐日积雪覆盖图像的去云装置的结构示意图。

具体实施方式

下面结合附图，对发明的具体实施方式作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

以下对本发明涉及的数据处理方法进行介绍：

1)基于积雪的时间连续性假设去云

此类方法主要依赖于积雪在短时间内不会快速消融的潜在假设，在时间尺度上实现去云。根据时间长短分为四种合成算法：上下午星数据合成、临近日分析、灵活数日结合和长时间序列分析。上下午星数据合成、临近日合成方法的使用有其局限性，前提假设是在一个连续的冰雪覆盖期或无雪期成立，而在积雪积累期和融化期可能每天都变。

2)基于积雪空间连续性假设去云

该类方法假设积雪或陆地在空间上连续而不破碎，根据云像素周围的非云像素信息将云重新分类，应用较多的是临近四像元法和临近八像元法。用临近四像元法去除了云，尽管去云量很少，但由于误差很小也构成去云方法的一部分，从而使所有步骤的总分类误差最小化。

3)基于高程的去云方法

该类方法基于这一事实：由于温度递减率和能量平衡，积雪在高海拔地区更早地积累，更晚地融化，积雪最低海拔之下没有积雪覆盖。此类方法包括三种典型算法：snowl方法、高程滤波和高程掩膜。其中Snowl方法的使用对云量有限制，当云量大于规定的阈值时，计算雪线和陆地线的高程会出现较大误差。阈值不同，方法的精度不同。当积雪空间异质性及破碎化比较严重时，snowl方法在山区的使用也受到一定的限制，黄晓东在青藏高原应用此方法时，首先进行了高程分带，以减小研究区高程差较大面积较广带来的影响，但积雪破碎化对算法产生的影响仍旧无法避免。

4)多传感器融合去云

多传感器融合法是利用被动微波不受云制约的优点，用AMSR-E雪水当量数据填补MODIS云区。然而，AMSR-E雪水当量产品对于积雪的判识具有较大误差，以美国西北部为测试区，对比每日AMSR-E数据和MODIS数据的分类精度发现，AMSR-E产品对积雪有着较大的低估误差，同时引入了更多的不确定因素。同年，在对阿拉斯加地区AMSR-E数据验证结果表明，数据精度只有68.5％，且存在对SWE的高估误差，在积雪积累期和融化期分类不太精确。在对中国北疆AMSR-E进行数据验证发现，在2002年11月至2005年3月的三个雪季中，总分类精度为65.1％～71.1％不等，积雪分类精度为65.6％～67.9％不等，AMSR-E雪水当量产品空间分辨率粗糙，尤其在积雪边缘地区低估了积雪面积。

5)其他数学方法

根据云周围像素的空间特征和地形特征来估计雪发生的概率。通过灵敏度和误差分析得到最小阈值，进而将其分类成积雪或陆地。Xia等以连续五天的MODIS积雪图像作为一个计算单元，用三维隐函数模拟这五天积雪边界的变化，可得到这段时间内每一时刻的积雪边界，界内即赋值为积雪，界外为无积雪覆盖的陆地。

图1示出了本发明一实施例提供的一种基于MODIS的逐日积雪覆盖图像的去云方法的流程示意图，包括：

S1、获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；

其中，所述第一预设时间段可根据具体的要求进行选择，例如，可选择365天；所述第一优选级可以为积雪>湖冰>湖泊>陆地。

在同一天中，高原地区云的出现和变化较快，而积雪覆盖情况不变，对MODIS数据进行上下午星积雪产品合成，参考Gao[17]的最大合成算法，公式表达如式(1)：

其中x,y分别表示像元S所在的行号和列号，t表示水文年的某一天，A和T分别表示下午星Aqua和上午星Terra，S(x,y,t)为像素值，合成优先级顺序定为：积雪>湖冰>湖泊>陆地，若合成图像仍有数据缺失或其他无效值，则赋值为云，与其他云像素一同进行下一步处理。对2009-2011两个积雪季每天的去云量进行统计，此步骤平均去除云量约11.0％，获得每日最大积雪合成产品。

S2、根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；

一般而言，阴天(多云)天气积雪接受太阳辐射较少，不易融化，降落在地面上会持续一段时间[31]，根据临近日的雪盖信息将当天的云像素重新分类。公式表达如式(2)：

S_(x,y,t)＝1if(S_(x,y,t-1)＝1and S_(x,y,t+1)＝1) (2)

其中x、y、t含义与式(1)中相同，若当日合成的积雪产品有云像素，前后日的相应位置同为积雪，则当日云像素分类为积雪，陆地覆盖同理。若前后日的相应位置至少一个为湖冰，则当日云像素赋值为湖冰，湖泊同理。其他情况则不去云。三天合成在2009-2011雪季中平均去除了7.5％的云量。

S3、将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于“永久性”陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；

S4、根据第二合成优先级，将连续八天的所述第三积雪覆盖合成图像进行“最大陆地”和“最大积雪”合成处理，得到八天合成处理结果，并根据所述八天处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；

基于以上高程滤波结果数据，开展8天最大积雪覆盖合成产品(类似MOIDS 8天积雪产品)和8天最大陆覆盖地合成，优先级顺序分别为：积雪>湖冰>湖泊>陆地、陆地>湖冰>湖泊>积雪，得到的合成结果分别记为图像a(8天最大积雪合成，也即8天最小陆地合成)、图像b(8天最大陆地合成，即8天最小积雪合成)。

积雪在8天之内逐渐融化或突然增加，在空间上有一定的渐变规律，每天的积雪面积均大于8天积雪的最小值积雪面积，而陆地面积大于8天陆地最小面积，算法用公式表示如式(4)、式(5)。图6显示了2011年1月9日-2011年1月17日这8天最小积雪面积与最小陆地面积，最小积雪和最小陆地面积之间存在的白色区域为不确定区域，即为8天积雪变化的范围。如此把每天的积雪动态范围以8天为时间窗口，控制在一定的动态变化范围，从而实现对于积雪去云过程的判识精度的控制。

其中，1和0分别表示积雪和陆地，表示8天最大积雪合成图像中第x行、第y列像素的值，类似。如果像素是云，任一最大合成图像中像素是湖冰或湖泊，则给云像素相应赋值为湖冰或湖泊。这一步骤去除云量最多，约16.5％。

S5、对所述第四积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

由于积雪覆盖与经度、纬度、坡度、坡向以及高程因素有关，且海拔高的地方比海拔低的地方被积雪覆盖的概率大。对整个区域用拟合预期雪线方法，算法描述如下：①提取雪线样本。积雪与陆地的分界线即为雪线，提取出雪线之上积雪的高程、坡度、坡向、经度、纬度。②拟合预期雪线高程。以高程为因变量，以其他四个影响因素为自变量，用多元线性回归拟合预期雪线高程。③将云像素重新分类。由云像素的坡度、坡向、经度、纬度可算出预期雪线高程，将其与云的实际高程比较，若实际高程较高，则将云赋值为积雪，反之赋值为陆地，公式表示如下：

其中，1和0分别表示积雪和陆地，为云像素对应的高程，表示计算出相同位置的预期雪线高程。由于青藏高原面积较大，且不同地区地形差异很大，按地形和高程将青藏高原分为7个子区域(如图6所示的算法流程图)，每个区分别使用拟合预期雪线方法去云。将得到结果与Google earth清晰的积雪图对比发现，验证显示，此方法基本将积雪区域和陆地区域正确的分离开来，拟合效果较好。这一步骤去除了最后剩余约5.1％的云覆盖，最后获得得到MODIS逐日无云积雪产品。

本实施例通过将多个卫星数据进行合成，并根据阈值对合成后的卫星图像进行多步去云操作，能够完全去除图像中影响动态监测的云，满足监测者的要求。

作为本实施例的可选方案，所述方法还包括：

S6、根据预设积雪覆盖精度指标，对所述MODIS无云积雪图像进行验证；

其中，所述预设积雪覆盖精度指标包括准确度、召回率、精度和f值。

为了评估和验证积雪的分类精度和效果，采用积雪面积精度指标来对每一步骤的积雪分类做验证评价，指标包括总分类精度(Accuracy)、召回率(Recall)、积雪分类精度(Precision)以及f值[33]。总精度即准确度，表示一个像素被正确分类的概率(可能性)，但是当图像中有大面积的无雪区域时，这一指标可能会带来误导，如若图像中只有10％积雪像元，那么总精度亦会由于大部分像元被分为无雪而达到90％。召回率表示检测到积雪像素的概率，即正确分类出的积雪像元数与实际所有积雪样本像素数的比值。积雪分类精度表示正确分类出的积雪像元数与分类出的所有积雪像元的比值。召回率和分类精度是相互影响的，召回率高表明正确分类的积雪像元数相对于地面站点观测的准确性较高，而精度高则表明正确分类的积雪像元数相对于MODIS中总积雪数的准确性较高，召回率高，并不代表分类精度高，因此采用f值(见公式13)或许是最合适的指标，这是召回率和积雪分类精度的调和平均值，能够较好的反映广泛无雪区域的缺失或错误分类的雪的情况。

具体地，所述基于“永久性”陆地和积雪覆盖的判别处理进一步包括：

S311、将所述第二积雪覆盖合成图像中，高程大于第一预设值的区域对应的像素点调整为积雪像素；

S312、若当前像素为云像素，所述云像素的高程大于第二预设值小于等于所述第一预设值，且所述云像素的云天数与积雪天数之和大于总天数与第三预设值的乘积，则将所述云像素调整为积雪像素；

S313、若所述云像素的云天数与陆地天数之和与所述总天数相同，且所述云天数小于所述总天数与第四预设值的乘积，则将所述云像素调整为陆地像素。

青藏高原山区分布着丰富的永久性冰雪，无论这些地方是否被云覆盖，都可直接将对应位置的云像素赋值为积雪，而有些地方在某些年份可能从不下雪或极少下雪，则可以通过判别而直接赋值为陆地。以2010-2011年积雪季为例，通过对比Google Earth和已发布的青藏高原永久性冰雪数据，并结合基于DEM的分布统计分析发现，高程在5800m以上的区域，可认为积雪覆盖(5800m以上区域面积包含于永久积雪面积)；由于积雪上面云层往往多于陆地，且避免MODIS积雪识别算法的错判可能，规定只要在一个雪季某像素积雪覆盖和云覆盖的时间综合大于全部时间的95％，即认为云可分类为积雪。类似地，如果某像素在这段时间内只出现陆地和云，但从未出现积雪，且云覆盖天数较少，则将云分类为陆地。经与Google Earth对照，云天数小于总天数的20％时，可保证精度不产生重大的损失。从而形成经验性的积雪去云的判断条件。

对每个时间段进行“永久”积雪和“永久”陆地的统计分析，判断条件如下：

①程>5800m的像素全部为“永久”积雪；

②高程在3000～5800之间的云像素，若满足：云天数+积雪天数>总天数*0.95，则云像素分类为积雪；

其中，为了避免出现云天数占大多数的情况，因此限定了云像素的高程在3000-5800，例如某像素一年中有300天为云像素，也符合云天数+积雪天数>总天数*0.95，但显示该云像素不一定为积雪像素；理论上，云天数+积雪天数＝总天数，但考虑到机器测量误差，因此采用云天数+积雪天数>总天数*0.95。

③若某像素满足：a.云天数+陆地天数＝总天数；b.云天数<总天数*0.2，则云像素分类成陆地。

具体地，所述临近四像元处理进一步包括：

S32、若所述云像素相邻的四个像素中至少三个为同类型非云像素，则将当前云像素调整为所述同类型非云像素；

其中，像素的数据类型包括：云、积雪、陆地、湖泊和湖冰；

所述同类型非云像素包括积雪像素，陆地像素、湖泊像素和湖冰像素。

青藏高原在积雪和融雪季节的积雪较为破碎，但对于多数情况而言，积雪在空间上连续的概率较大，可根据云像素周围的非云像素信息对其进行分类。如果与云像素相邻最近的四个像元中，至少三个为积雪，那么将中心云像素赋值为雪；如果至少三个为陆地，那么将中心云像素赋值为陆地；湖泊、湖冰同理。该方法去除云量较少，为5％，但分类精度很高，为整体去云算法的使用减少累积误差。

具体地，所述高程滤波处理进一步包括：

S33、若当前像素为积雪像素，所述积雪像素相邻的八个像素中有云像素，且所述云像素对应的高程大于所述积雪像素的高程，则将所述云像素调整为积雪像素。

一般而言，海拔较高地方温度较低，积雪容易积累，根据这一原理使用高程滤波方法，利用云像素周围信息及SRTM DEM辅助数据对上一步完成的MODIS逐日积雪图像进行去云。公式描述如式(3)：

S_(x+i,y+j,t)＝1if S_(x,y,t)＝1and H_(x+i,y+j)＞H_(x,y)，i、j＝±1 (3)

以积雪像素为中心，若与其相邻的八个像元中有云，且云像素高程大于中心积雪高程，则将云像素赋值为积雪，该步骤去除0.9％的云量，精度可以获得保证。

具体地，所述将高山冰川覆盖区的湖水误判为阴影区处理进一步包括：

S34、根据现有湖泊辅助数据像素的矢量边界，将所述矢量边界之外的湖泊像素和湖冰像素均调整为云像素。

通过对比世界湖泊矢量边界，发现MODIS数据在青藏高原的西部地区，常把阴影区错误地判断为湖冰或湖泊，对后面的去云过程造成误差累计。为减小这种误差，用已有湖泊矢量边界作为辅助数据，对错判的湖泊湖冰进行修正。边界线之内的湖泊或湖冰保留，界线之外错分的湖泊和湖冰重新赋值为云，进行下一步判断。

具体地，所述拟合预期雪线方法进一步包括：

S51、根据云像素的坡度、坡向、经度和纬度，得到云像素的高程曲线，并根据云的实际高程，将小于实际高程的云像素调整为积雪像素。

由于积雪覆盖与经度、纬度、坡度、坡向以及高程因素有关，且海拔高的地方比海拔低的地方被积雪覆盖的概率大。对整个区域用拟合预期雪线方法，算法描述如下：①提取雪线样本。积雪与陆地的分界线即为雪线，提取出雪线之上积雪的高程、坡度、坡向、经度、纬度。②拟合预期雪线高程。以高程为因变量，以其他四个影响因素为自变量，用多元线性回归拟合预期雪线高程。③将云像素重新分类。由云像素的坡度、坡向、经度、纬度可算出预期雪线高程，将其与云的实际高程比较，若实际高程较高，则将云赋值为积雪，反之赋值为陆地。亚洲地区的高山区(如青藏高原)地形复杂，引起热量水汽分布极其不均，因此导致积雪的时空分布极不均匀。青藏高原地区独特复杂的地形以及积雪破碎的情况，使得当前的每日无云积雪遥感动态监测和估算存在不确定性，需要探索不同的算法。

目前已有的方法中往往在保证精度的条件下，无法获得逐日无云积雪产品，如可变天数产品，采用上午星数据MOD10A1测试积雪2～5天的连续性、三天合成、五天合成等，只要积雪产品达到云量小于10％则停止工作。而对于获得逐日无云积雪产品，其他各种方法很难适用到青藏高原的复杂地形区域，而多传感器融合方法，选择某天AMSR-E雪水当量产品和MODIS每日积雪产品，忽略了AMSR-E对于积雪面积判别的误差。

本实施例通过研究以MODIS逐日积雪产品数据(MOD10A1和MYD10A1)为基础，综合考虑青藏高原的地形特点和积雪分布特征，在通过针对不同的算法步骤的基础上，逐步控制剩余云量和精度，形成新的去云算法流程，并最终得到MODIS逐日无云积雪产品系列算法的开发，为青藏高原积雪动态监测和历史积雪参数数据的重建提供了参考。

青藏高原地区面积广阔，复杂多样的地形及积雪空间异质性大的特点决定了去云方法的研究需考虑区域的适用性，根据对青藏高原的积雪分布特征的分区假设的检验研究，发现青藏高原地区云覆盖像元的再分类需要考虑，1)高原腹地，由于对流过程的复杂和地形坡度的较缓，积雪的出现时间短，且变化较快。2)高原周边山区的积雪存在永久性冰雪问题，且积雪受高程的控制明显(温度较低)，降水的概率受高原西风和季风区的影响原因不同，积雪的分布特征不同。因此，有关数字高程模型的算法在青藏高原应采用分区的方式，在腹地其效果将产生较大的误差。

本实施例采用步骤A1-A8对青藏高原MODIS积雪产品进行去云，逐步减少云量，最终得到逐日无云产品。算法流程图见图2，其过程描述和分析如下：

A1、MODIS上下午星数据合成

在同一天中，高原地区云的出现和变化较快，而积雪覆盖情况不变，对MODIS数据进行上下午星积雪产品合成，参考Gao[17]的最大合成算法，公式表达如式(1)：

其中x,y分别表示像元S所在的行号和列号，t表示水文年的某一天，A和T分别表示下午星Aqua和上午星Terra，S(x,y,t)为像素值，合成优先级顺序定为：积雪>湖冰>湖泊>陆地，若合成图像仍有数据缺失或其他无效值，则赋值为云，与其他云像素一同进行下一步处理。对2009-2011两个积雪季每天的去云量进行统计，此步骤平均去除云量约11.0％，获得每日最大积雪合成产品。

A2、三天合成

一般而言，阴天(多云)天气积雪接受太阳辐射较少，不易融化，降落在地面上会持续一段时间[31]，根据临近日的雪盖信息将当天的云像素重新分类。公式表达如式(2)：

S_(x,y,t)＝1if(S_(x,y,t-1)＝1and S_(x,y,t+1)＝1) (2)

其中x、y、t含义与式(1)中相同，若当日合成的积雪产品有云像素，前后日的相应位置同为积雪，则当日云像素分类为积雪，陆地覆盖同理。若前后日的相应位置至少一个为湖冰，则当日云像素赋值为湖冰，湖泊同理。其他情况则不去云。三天合成在2009-2011雪季中平均去除了7.5％的云量。

A3、“永久积雪、陆地”判别法

青藏高原山区分布着丰富的永久性冰雪，无论这些地方是否被云覆盖，都可直接将对应位置的云像素赋值为积雪，而有些地方在某些年份可能从不下雪或极少下雪，则可以通过判别而直接赋值为陆地。以2010-2011年积雪季为例，通过对比Google Earth和已发布的青藏高原永久性冰雪数据，并结合基于DEM的分布统计分析发现，高程在5800m以上的区域，可认为积雪覆盖(5800m以上区域面积包含于永久积雪面积)；由于积雪上面云层往往多于陆地，且避免MODIS积雪识别算法的错判可能，规定只要在一个雪季某像素积雪覆盖和云覆盖的时间综合大于全部时间的95％，即认为云可分类为积雪。类似地，如果某像素在这段时间内只出现陆地和云，但从未出现积雪，且云覆盖天数较少，则将云分类为陆地。经与Google Earth对照，云天数小于总天数的20％时，可保证精度不产生重大的损失。从而形成经验性的积雪去云的判断条件。

对每个时间段进行“永久”积雪和“永久”陆地的统计分析，判断条件如下：

①高程>5800m的像素全部为“永久”积雪；

②高程在3000～5800之间的云像素，若满足：云天数+积雪天数>总天数*0.95，则云像素分类为积雪；

③若某像素满足：a.云天数+陆地天数＝总天数；b.云天数<总天数*0.2，则云像素分类成陆地。

如图3为MODIS逐日无云积雪产品永久积雪、陆地判别算法流程图。如图4所示，用此图像中的“永久”积雪和“永久”陆地代替这一雪季每日对应位置的云像素，统计此步骤在2009-2011两个雪季的去云量，平均值达3.1％。

A4、临近四像元法

如图5示出了本实施例提供的MODIS逐日无云积雪产品临近四像元算法流程图。青藏高原在积雪和融雪季节的积雪较为破碎，但对于多数情况而言，积雪在空间上连续的概率较大，可根据云像素周围的非云像素信息对其进行分类。如果与云像素相邻最近的四个像元中，至少三个为积雪，那么将中心云像素赋值为雪；如果至少三个为陆地，那么将中心云像素赋值为陆地；湖泊、湖冰同理。该方法去除云量较少，为5％，但分类精度很高，为整体去云算法的使用减少累积误差。

A5、高程滤波法

一般而言，海拔较高地方温度较低，积雪容易积累，根据这一原理使用高程滤波方法，利用云像素周围信息及SRTM DEM辅助数据对上一步完成的MODIS逐日积雪图像进行去云。公式描述如式(3)：

S_(x+i,y+j,t)＝1if S_(x,y,t)＝1and H_(x+i,y+j)＞H_(x,y)，i、j＝±1 (3)

以积雪像素为中心，若与其相邻的八个像元中有云，且云像素高程大于中心积雪高程，则将云像素赋值为积雪，该步骤去除0.9％的云量，精度可以获得保证。

A6、修改阴影区分类错误

通过对比世界湖泊矢量边界，发现MODIS数据在青藏高原的西部地区，常把阴影区错误地判断为湖冰或湖泊，对后面的去云过程造成误差累计。为减小这种误差，用已有湖泊矢量边界作为辅助数据，对错判的湖泊湖冰进行修正。边界线之内的湖泊或湖冰保留，界线之外错分的湖泊和湖冰重新赋值为云，进行下一步判断。

A7、八天最大积雪、最大陆地面积掩膜法

基于以上高程滤波结果数据，开展8天最大积雪覆盖合成产品(类似MOIDS 8天积雪产品)和8天最大陆覆盖地合成，优先级顺序分别为：积雪>湖冰>湖泊>陆地、陆地>湖冰>湖泊>积雪，得到的合成结果分别记为图像a(8天最大积雪合成，也即8天最小陆地合成)、图像b(8天最大陆地合成，即8天最小积雪合成)。

积雪在8天之内逐渐融化或突然增加，在空间上有一定的渐变规律，每天的积雪面积均大于8天积雪的最小值积雪面积，而陆地面积大于8天陆地最小面积，算法用公式表示如式(4)、式(5)。图6显示了2011年1月9日-2011年1月17日这8天最小积雪面积与最小陆地面积，最小积雪和最小陆地面积之间存在的白色区域为不确定区域，即为8天积雪变化的范围。如此把每天的积雪动态范围以8天为时间窗口，控制在一定的动态变化范围，从而实现对于积雪去云过程的判识精度的控制。

其中，1和0分别表示积雪和陆地，表示8天最大积雪合成图像中第x行、第y列像素的值，类似。如果像素是云，任一最大合成图像中像素是湖冰或湖泊，则给云像素相应赋值为湖冰或湖泊。这一步骤去除云量最多，约16.5％。

A8、拟合预期雪线方法

由于积雪覆盖与经度、纬度、坡度、坡向以及高程因素有关，且海拔高的地方比海拔低的地方被积雪覆盖的概率大。对整个区域用拟合预期雪线方法，算法描述如下：①提取雪线样本。积雪与陆地的分界线即为雪线，提取出雪线之上积雪的高程、坡度、坡向、经度、纬度。②拟合预期雪线高程。以高程为因变量，以其他四个影响因素为自变量，用多元线性回归拟合预期雪线高程。③将云像素重新分类。由云像素的坡度、坡向、经度、纬度可算出预期雪线高程，将其与云的实际高程比较，若实际高程较高，则将云赋值为积雪，反之赋值为陆地，公式表示如下：

其中，1和0分别表示积雪和陆地，为云像素对应的高程，表示计算出相同位置的预期雪线高程。由于青藏高原面积较大，且不同地区地形差异很大，按地形和高程将青藏高原分为7个子区域(如图6所示的算法流程图)，每个区分别使用拟合预期雪线方法去云。将得到结果与Google earth清晰的积雪图对比发现，验证显示，此方法基本将积雪区域和陆地区域正确的分离开来，拟合效果较好。这一步骤去除了最后剩余约5.1％的云覆盖，最后获得得到MODIS逐日无云积雪产品。

如图7为本实施例提供的MODIS逐日无云积雪产品拟合雪线算法流程图。

经过以上8个步骤，逐步实现了整个去云算法的过程，形成基于经验参数的青藏高原逐日无云积雪产品算法流程。

测试结果与分析

B1、各方法步骤去云效果分析

分析研究区内2009-2011两个积雪季，共计424个时相，在整个去云过程中云覆盖的动态变化，如图8所示，MYD10A1数据和MOD10A1云量较大，平均分别达到47.1％和41.7％，云量超过50％的分别有186天和121天，且多分布在每年的12月至次年3月，10月、11月和次年4月、5月云量相对减少，云的大量存在严重制约着积雪动态监测。上下午星合成和三天连续合成分别去除了11.0％和7.5％的云覆盖，去云量较多。“永久积雪、陆地”判别方法涉及到长时间的积雪覆盖信息，去除了永久积雪和永久陆地上覆盖的云，约为3.1％。临近四像元法去除零星分散约0.5％的云像素，去云量较少但精度较高。高程滤波方法将积雪边缘海拔较高的云像素重新分类，去除约0.9％的云量。8天最大积雪、最大陆地范围掩膜方法利用短时间内雪盖信息的相关性，将16.4％的云像素重新分类为积雪或陆地像素，去云量最多，但是其内在的关系是连续多天的积雪的不变性假设为基础，是三天连续积雪去云方法的一种延续。通过以上七个步骤，MODIS积雪合成产品仅剩5.1％的云覆盖，采用拟合预期雪线算法去除剩余全部的云，或者逐日无云积雪产品。图9示出了本实施例提供的2009-2011年积雪季过程产品云量分析图。

B2、精度验证

为了评估和验证积雪的分类精度和效果，采用积雪面积精度指标来对每一步骤的积雪分类做验证评价，指标包括总分类精度(Accuracy)、召回率(Recall)、积雪分类精度(Precision)以及f值[33]。总精度即准确度，表示一个像素被正确分类的概率(可能性)，但是当图像中有大面积的无雪区域时，这一指标可能会带来误导，如若图像中只有10％积雪像元，那么总精度亦会由于大部分像元被分为无雪而达到90％。召回率表示检测到积雪像素的概率，即正确分类出的积雪像元数与实际所有积雪样本像素数的比值。积雪分类精度表示正确分类出的积雪像元数与分类出的所有积雪像元的比值。召回率和分类精度是相互影响的，召回率高表明正确分类的积雪像元数相对于地面站点观测的准确性较高，而精度高则表明正确分类的积雪像元数相对于MODIS中总积雪数的准确性较高，召回率高，并不代表分类精度高，因此采用f值(见公式13)或许是最合适的指标，这是召回率和积雪分类精度的调和平均值，能够较好的反映广泛无雪区域的缺失或错误分类的雪的情况。

青藏高原2009年10月1日-2011年4月30日两个雪季共424个时相，145个地面台站提供每日雪深数据，共61480对估算和观测值。通过混淆矩阵及式(10)、(11)、(12)、(13)对积雪分类图像进行评价，分析晴空天气下MODIS数据本身(MOD10A1、MYD10A1)以及去云过程中各合成产品的精度，在精度评价中考虑以下几个方面的样本数：1)地面台站记录和积雪图像均有雪(雪深>0)的样本数a；2)地面台站记录有雪而图像分为无雪的样本数b，即漏测数；3)地面台站记录无雪而图像分为有雪类型的样本数c；4)地面台站记录和图像均无雪的样本数d。

表1精度验证混淆矩阵

积雪精度指标公式如下：

原始数据MOD10A1、MYD10A1及各过程产品在混淆矩阵中的样本数统计如表2，在积雪去云的处理过程中，修改湖泊、湖冰错判改变的像元数少，故省略统计结果。青藏高原积雪大多较薄，且分布斑驳分散，高原地区下垫面异质性较大，地面台站雪深数据与MODIS积雪图像像元之间空间尺度的差异会造成精度验证的误差。当融雪期积雪较薄(如小于3cm)时，可能出现这种情况，包括地面台站所在的小面积区域观测为有雪，而MODIS图像500m×500m像素判断为无雪，从而产生积雪的低估误差。也可能台站处于城市，由于温度高于四周而较快融雪，此时MODIS图像判断为有雪，如此出现积雪高估误差。另外，对于斑驳积雪的深度较小的情况，光学可以做一定的穿透，且混合像元问题较为严重，因此采用一定深度的积雪开展验证，将更好地对积雪去云的判别过程开展评价。若只保留积雪深度大于3cm的样本数据，图像和地面观测同时记录为积雪的样本数记为a3，图像记录为积雪而地面记录为陆地的样本数记录为b3，仍用上式计算，各合成产品的分类精度较采用全部积雪观测数据有所提高。

MODIS数据本身的总分类精度在98％以上，当积雪深度大于3cm时，其f值超过90％，在台站所观测的区域精度较高。从整个积雪去云的角度来看，积雪的召回率和积雪分类精度直到少云步骤(步骤七)的过程中，在两个冬季期间的积雪判别中，其精度基本上没有损失，主要的原因一方面是由于积雪去除的过程中判别的精度损失较少，也存在着部分在在有积雪观测的区域，并未能受到太多的影响，而对于最不稳定的去云的最后一步采用拟合线的时候获得了较大的误差，但是依然可达86.31％的f值。

总体来看，由表2可见，高程滤波之前的去云步骤对分类精度影响不大，上下午最大积雪合成和高程滤波方法稍微提高了积雪分类精度，最大积雪、陆地面积掩膜方法使合成产品的精度有所将低，但去云量最多。拟合预期雪线方法精度最低，用于去云最后一步，但当雪深大于3cm时，无云积雪产品的总分类精度、召回率、积雪分类精度及f值仍可达96.6％、89％、83.8％和86.3％，因此积雪分类上依然存在一定的提升空间，但是这种提高需要考虑本身MODIS积雪产品在高原地区地表不均一，积雪斑驳和动态变化快等引起的不足。

表2各算法过程产品精度一览表

青藏高原积雪平均于每年9月份开始，次年4-5月份结束。地形复杂多样，积雪分布不均，海拔较高的山区积雪开始较早，结束较晚，积雪持续时间长。永久积雪(冰川)的面积程减小趋势，转化为季节性积雪，永久陆地面积程增多趋势，这与全球变暖有关，且对气候起着反馈与调节作用。

通过多次测试并与Google earth较清晰的图像对比，得到如下“永久”积雪和“永久”陆地判别条件：①高程>5800m的像素全部为“永久”积雪；②高程在3000～5800之间的云像素，若满足：云天数+积雪天数>总天数*0.95，则该像素为“永久”积雪；③若某像素满足：a.云天数<总天数*0.2，b.云天数+陆地天数＝总天数，则云像素分类成陆地。如果与云像素相邻最近的四个像元中，至少三个为积雪，那么将中心云像素赋值为雪；如果至少三个为陆地，那么将中心云像素赋值为陆地；湖泊、湖冰同理。以积雪像素为中心，若与其相邻的八个像元中有云，且云像素高程大于中心积雪高程，则将云像素赋值为积雪。湖泊和湖冰为连续8天的最大值，认为最大积雪图像中的陆地面积是8天最小，最大陆地图像中的积雪面积是8天最小，两个图像融合得到最小积雪面积、最小陆地面积图像。考虑这8天积雪覆盖在时间上相关，在空间上有一定的渐变规律，其中任意一天的积雪和陆地面积均大于8天的最小值，可用这一最小面积代替每天相应位置的云像素。

按照去云原理，对常用的去云算法进行分类，分析每一类方法在青藏高原地区的适用性。

本实施例研究区为青藏高原，将常用的去云方法按原理分成5大类，验证每种方法本身的正确率，进而确定是否可用于青藏高原。分析每一种去云算法的误差来源，在青藏高原适用性较好的方法继续使用，对于可改进的算法则根据青藏高原的具体情况予以改进，然后应用。结合青藏高原特殊的地形地貌，提出适用于青藏高原独特的去云方法，以便去除更多的云，并尽可能降低误差累积，最后得到完整的去云算法。

青藏高原的逐日无云积雪的新型算法流程，算法流程实现的八个步骤的有机整合。

通过本实施例不但能够加快青藏高原逐日无云积雪的计算，为后来研究者的提供了算法途径，可以满足当前研究的需求。本发明的研究表明在高原地区，当积雪深度＞3cm时，无云积雪产品总分类精度达到96.6％，积雪分类精度达89.0％，整个算法流程对MODIS积雪产品去云的精度损失较低，数据可靠性较高。

青藏高原地形复杂多样，积雪分布不均，海拔较高的山区积雪分布较多，而高原腹地降雪很少。根据积雪的时间连续性、空间连续性，以及雪盖与高程的关系，采用不同方法进行去云，获得了精度较好的MODIS无云积雪产品，并可供一定的模型作为数据，当前已形成可供使用的产品数据集。然而，针对青藏高原的山区的积雪面积的监测精度的进一步提高，还需要基于原始的观测值的基础之上开展非均匀复杂地形条件下混合像元的分析，并充分考虑高山区的观测和入射角等因素，从根本上提高MODIS数据的积雪分类精度。台站分布不均及其代表性也影响着无云产品的精度验证，可考虑采用更多不同的高分辨率数据的方式，开展积雪精度产品的验证，以期获得实时监测青藏高原积雪动态变化、防灾减灾、气候预测等应用。

通过不断修改青藏高原云效果的参数，结合青藏高原地形和山地积雪特征，利用8个步骤去除了所有云覆盖，并选取2009-2011年度两个积雪季为研究期，采用该期间176个地面台站观测数据，对各合成步骤及最后结果进行精度验证。结果表明，MODIS每日无云积雪产品总分类精度较高，为94.7％。考虑到光学遥感对薄雪的探测精度本身较差，若去掉雪深<＝3cm的样本，总分类精度达96.6％。算法对MODIS积雪产品去云效果较好，用此方法批量生产青藏高原2002-2015年的无云积雪产品。

MYD10A1数据和MOD10A1云量较大，平均分别达到47.1％和41.7％，通过上下午星合成和三天连续合成分别去除了11.0％和7.5％的云覆盖，去云量较多。“永久积雪、陆地”判别方法涉及到长时间的积雪覆盖信息，去除了永久积雪和永久陆地上覆盖的云，约为3.1％。临近四像元法去除零星分散约0.5％的云像素，去云量较少但精度较高。高程滤波方法将积雪边缘海拔较高的云像素重新分类，去除约0.9％的云量。8天最大积雪、最大陆地范围掩膜方法利用短时间内雪盖信息的相关性，将16.4％的云像素重新分类为积雪或陆地像素，去云量最多，但是其内在的关系是连续多天的积雪的不变性假设为基础，是三天连续积雪去云方法的一种延续。通过以上七个步骤，MODIS积雪合成产品仅剩5.1％的云覆盖，采用拟合预期雪线算法去除剩余全部的云，或者逐日无云积雪产品。

图10示出了本实施例提出的一种基于MODIS的逐日积雪覆盖图像的去云装置的结构示意图，所述装置包括上下午合成模块11、三天合成模块12、处理后合成模块13、八天合成模块14和雪线拟合模块15，其中，

上下午合成模块11用于获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；

三天合成模块12用于根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；

处理后合成模块13用于将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于“永久性”陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；

八天合成模块14用于根据第二合成优先级，将连续八天的所述第三积雪覆盖合成图像进行“最大陆地”和“最大积雪”合成处理，得到八天合成处理结果，并根据所述八天处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；

雪线拟合模块15用于对所述第三积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

具体地，上下午合成模块11获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；三天合成模块12根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；处理后合成模块13将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于“永久性”陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；八天合成模块14根据第三合成优先级，将连续八天的所述第三积雪覆盖合成图像进行“最大陆地”和“最大积雪”合成处理，得到八天合成处理结果，并根据所述八天处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；雪线拟合模块15对所述第四积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

本实施例通过将多个卫星数据进行合成，并根据阈值对合成后的卫星图像进行多步去云操作，能够完全去除图像中影响动态监测的云，满足监测者的要求。

进一步地，所述装置还包括：

图像验证模块16，用于根据预设积雪覆盖精度指标，对所述MODIS无云积雪图像进行验证；

其中，所述预设积雪覆盖精度指标包括准确度、召回率、精度和f值。

更进一步地，所述处理后合成模块中积雪覆盖的判别处理包括：

第一调整单元，用于将所述第二积雪覆盖合成图像中，高程大于第一预设值的区域对应的像素点调整为积雪像素；

第二调整单元，用于若当前像素为云像素，所述云像素的高程大于第二预设值小于等于所述第一预设值，且所述云像素的云天数与积雪天数之和大于总天数与第三预设值的乘积，则将所述云像素调整为积雪像素；

第三调整单元，用于若所述云像素的云天数与陆地天数之和与所述总天数相同，且所述云天数小于所述总天数与第四预设值的乘积，则将所述云像素调整为陆地像素。

本实施例所述的基于MODIS的逐日积雪覆盖图像的去云装置可以用于执行上述方法实施例，其原理和技术效果类似，此处不再赘述。

本发明的说明书中，说明了大量具体细节。然而，能够理解，本发明的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。

Claims (10)

1.一种基于MODIS的逐日积雪覆盖图像的去云方法，其特征在于，包括：

获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；

根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；

将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于永久性陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；

根据第二合成优先级，将连续八天的所述第三积雪覆盖合成图像进行最大陆地和最大积雪合成处理，得到八天合成处理结果，并根据所述八天合成处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；

对所述第四积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

2.根据权利要求1所述的方法，其特征在于，还包括：

根据预设积雪覆盖精度指标，对所述MODIS无云积雪图像进行验证；

其中，所述预设积雪覆盖精度指标包括准确度、召回率、精度和f值。

3.根据权利要求1所述的方法，其特征在于，所述基于永久性陆地和积雪覆盖的判别处理进一步包括：

将所述第二积雪覆盖合成图像中，高程大于第一预设值的区域对应的像素点调整为积雪像素；

若当前像素为云像素，所述云像素的高程大于第二预设值小于等于所述第一预设值，且所述云像素的云天数与积雪天数之和大于总天数与第三预设值的乘积，则将所述云像素调整为积雪像素；

若所述云像素的云天数与陆地天数之和与所述总天数相同，且所述云天数小于所述总天数与第四预设值的乘积，则将所述云像素调整为陆地像素。

4.根据权利要求3所述的方法，其特征在于，所述临近四像元处理进一步包括：

若所述云像素相邻的四个像素中至少三个为同类型非云像素，则将当前云像素调整为所述同类型非云像素；

其中，像素的数据类型包括：云、积雪、陆地、湖泊和湖冰；

所述同类型非云像素包括积雪像素，陆地像素、湖泊像素和湖冰像素。

5.根据权利要求4所述的方法，其特征在于，所述高程滤波处理进一步包括：

若当前像素为积雪像素，所述积雪像素相邻的八个像素中有云像素，且所述云像素对应的高程大于所述积雪像素的高程，则将所述云像素调整为积雪像素。

6.根据权利要求5所述的方法，其特征在于，所述将高山冰川覆盖区的湖水误判为阴影区处理进一步包括：

根据现有湖泊辅助数据像素的矢量边界，将所述矢量边界之外的湖泊像素和湖冰像素均调整为云像素。

7.根据权利要求1所述的方法，其特征在于，所述拟合预期雪线方法进一步包括：

根据云像素的坡度、坡向、经度和纬度，得到云像素的高程曲线，并根据云的实际高程，将小于实际高程的云像素调整为积雪像素。

8.一种基于MODIS的逐日积雪覆盖图像的去云装置，其特征在于，包括：

上下午合成模块，用于获取第一预设时间段内MODIS的上午星和下午星逐日积雪覆盖图像，并根据第一合成优先级将每天的所述上午星和所述下午星的积雪覆盖图像进行合成，得到第一积雪覆盖合成图像；

三天合成模块，用于根据第一合成规则将连续三天的所述第一积雪覆盖合成图像进行合成，得到第二积雪覆盖合成图像；

处理后合成模块，用于将所述第一预设时间段划分为若干个第二预设时间段，根据第二合成规则对所述第二预设时间段内的所述第二积雪覆盖合成图像分别进行基于永久性陆地和积雪覆盖的判别处理、临近四像元处理、高程滤波处理和将高山冰川覆盖区的湖水误判为阴影区处理，得到第三积雪覆盖合成图像；

八天合成模块，用于根据第二合成优先级，将连续八天的所述第三积雪覆盖合成图像进行最大陆地和最大积雪合成处理，得到八天合成处理结果，并根据所述八天合成处理结果对所述第三积雪覆盖合成图像进行云覆盖像素的积雪判别，得到第四积雪覆盖合成图像；

雪线拟合模块，用于对所述第四积雪覆盖合成图像进行区域划分，根据拟合预期雪线方法对每个划分后的区域进行去云处理，得到完全去云的MODIS无云积雪图像。

9.根据权利要求8所述的装置，其特征在于，还包括：

图像验证模块，用于根据预设积雪覆盖精度指标，对所述MODIS无云积雪图像进行验证；

其中，所述预设积雪覆盖精度指标包括准确度、召回率、精度和f值。

10.根据权利要求8所述的装置，其特征在于，所述处理后合成模块中基于永久性陆地和积雪覆盖的判别处理进一步包括：

第一调整单元，用于将所述第二积雪覆盖合成图像中，高程大于第一预设值的区域对应的像素点调整为积雪像素；

第二调整单元，用于若当前像素为云像素，所述云像素的高程大于第二预设值小于等于所述第一预设值，且所述云像素的云天数与积雪天数之和大于总天数与第三预设值的乘积，则将所述云像素调整为积雪像素；

第三调整单元，用于若所述云像素的云天数与陆地天数之和与所述总天数相同，且所述云天数小于所述总天数与第四预设值的乘积，则将所述云像素调整为陆地像素。