CN104657952A - 一种结合时空信息的遥感逐日积雪产品去云方法 - Google Patents

Abstract

本发明公开了一种结合时空信息的遥感逐日积雪产品去云方法，首先根据DEM数据和土地利用数据对遥感逐日积雪产品进行积雪单元划分；再利用DEM数据和多时相数据分别得到每个云像元在空间和时间上的积雪概率；最后利用加权平均的方法获得云像元的积雪概率，评判云像元是积雪还是无雪。本发明充分结合了积雪在空间上和时间上的相关性，以及地形因子对积雪的影响，有效地减少云对积雪产品的影响，提高了积雪产品的精度，并且算法简单，计算效率高，易于实现。

Description

一种结合时空信息的遥感逐日积雪产品去云方法

技术领域

本发明属于遥感图像处理技术领域，涉及一种结合时空信息的遥感逐日积雪产品去云方法，利用积雪的特性及地形信息，有效地恢复积雪产品中云像元信息。

背景技术

由于遥感逐日积雪产品有较高的时间分辨率和空间分辨率，且在晴空条件下与地面实测数据吻合精度高，因此，遥感逐日积雪产品被广泛应用于积雪的检测中。然而，云量高是限制其应用的一个主要因素。目前，常见的去云方法有：两卫星合成，相邻时间推导，空间滤波，基于雪线的方法，根据积雪过渡期确定阈值的方法，多天合成(包括固定天数和非固定天数合成)等，采用其中一种方法或多种方法依次去云。其中，两卫星合成和相邻时间推导的方法运用比较广泛，去云量较大，精度高，但不能完全解决云覆盖问题；空间滤波方法主要基于4邻域或8邻域的信息去云，去云量小，且牺牲了部分空间分辨率；基于雪线的方法，随着云覆盖量增大，精度降低；根据积雪过渡期确定阈值的方法，只考虑了一个积雪周期，忽略了积雪累积和融化期的交替变化；多天合成方法，包括固定天数和可变天数合成，虽然能有效减少云覆盖量，但降低了时间分辨率。

总之，上述去云方法仅考虑到时间或空间单方面对积雪的影响，故在某些情况下不能完全去除云覆盖，并且对遥感逐日积雪产品的修复效果不佳。

发明内容

针对现有技术的上述不足，本发明提供了一种结合时空信息的遥感逐日积雪产品去云方法，充分利用积雪在空间和时间上的特性以及地形因子对积雪的影响，采用加权平均方法计算积雪概率判断云像元的类别。

本发明所采用的技术方案是：一种结合时空信息的遥感逐日积雪产品去云方法，其特征在于，包括以下步骤：

步骤1：原始的遥感逐日积雪产品重分类为积雪、无雪和云三类，同时根据DEM数据和土地利用数据对遥感逐日积雪产品进行积雪单元划分；

步骤2：利用DEM数据和多时相数据分别计算每个云像元在空间和时间上的积雪概率；

步骤3：利用加权平均方法获得云像元的积雪概率，评判云像元是积雪还是无雪。

作为优选，步骤1中所述的对遥感逐日积雪产品进行积雪单元划分，划分后将云量大于70％的积雪单元进行合并，形成新的积雪单元。

作为优选，步骤2中所述的计算每个云像元在空间和时间上的积雪概率，具体实现过程是：

首选在每个积雪单元的高程中，找到雪线H_S和无雪线H_L；其中，高于雪线H_S的无云像元都为积雪，低于无雪线H_L的无云像元都为无雪；在H_S、H_L之间等间距划分高程，并且高程大于等于H_S和小于等于H_L各为一类；

然后，在每个高程带中，根据南、北和东西三类坡向继续划分高程；

在高程、坡向分类基础上，每个云像元在其所在类中的积雪概率为P_H，具体公式如下：

P_H＝N_s/N_nocloud    (1)；

其中，N_s是该类中积雪像元数目，N_nocloud是该类中非云像元数目；

P_T是当前云像元在时间上前后各5个非云像元中的积雪概率，公式如下：

$P_T=\frac{\mathrm{\Σ}1/d_{\mathrm{index}}}{\underset{i=1}{\overset{10}{\mathrm{\Σ}}}1/d_i},\mathrm{index}\∈\{\mathrm{1,2,3},...,10\}---\left(2\right);$

其中，d_i是非云像元i到当前云像元的时相距离，d_index是非云像元中积雪像元到当前云像元的时相距离。

作为优选，步骤3中的具体实现过程是，根据下列公式得到云像元积雪的总体概率P_s，当P_s大于等于P₀时，将该云像元判断为积雪，反之为无雪：

P_s＝w₁×P_H+w₂×P_T       (3)；

其中，P₀为设定的阈值，w₁和w₂分别是P_H和P_T的权重，且w₁+w₂＝1。

作为优选，步骤3中所述的w₁和w₂，其计算方法是依据积雪单元中的非云像元的信息获取，具体步骤如下：

首先，在每个积雪单元中，利用步骤3中的方法求得非云像元在空间和时间上的积雪概率；

然后，在[0,1]区间内，以0.01为步长设定w₁，以公式(3)求得非云像元的积雪概率，判断其为积雪还是无雪，并将判断结果视为非云像元的预测值；

最后，将非云像元的原始值作为真实值，算出所述的遥感逐日积雪产品去云方法在非云像元上的精度oa，公式如下所示：

$\mathrm{oa}=\frac{N_s^s+N_{\mathrm{ns}}^{\mathrm{ns}}}{N_{\mathrm{nc}}}---\left(4\right);$

其中，为积雪单元中非云像元的预测值和真实值都为积雪的数目，为积雪单元中非云像元的预测值和真实值都为无雪的数目，N_nc为积雪单元中非云像元的总数目。

本发明的特色在于，能很好的结合积雪在空间上与流域、土地利用等具有相关性的特点和在时间上持续性的特点，并加入高程和坡向对积雪的影响，如阳坡的融雪速度快于阴坡；随着高程增大，积雪覆盖面积增大等。同时，本发明能生成完全无云的高精度积雪产品，且保持积雪产品的空间分辨率和时间分辨率不变，对积雪变化研究有着重要应用价值。

附图说明

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

具体实施方式

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

请见图1，本实施例提供的一种结合时空信息的遥感逐日积雪产品去云方法，包括以下步骤：

步骤1：将原始的遥感逐日积雪产品重分类为积雪、无雪和云三类，同时根据DEM数据和土地利用数据对遥感逐日积雪产品进行积雪单元划分；

首先，利用DEM数据和水系图将遥感影像所在区域划分为多个子流域。然后在子流域的基础上，根据土地利用数据将区域划分为多个积雪单元。为了保持每个积雪单元的云量小于70％，将云量大于70％的积雪单元与周围的积雪单元合并，直至新形成的积雪单元的云量小于70％。若最后所有的积雪单元都合并在一起了，则将整个区域视为一个积雪单元。

步骤2：利用DEM数据和多时相数据分别计算每个云像元在空间和时间上的积雪概率；具体实现过程是：

首选在每个积雪单元的高程中，找到雪线H_S和无雪线H_L；其中，高于雪线H_S的无云像元都为积雪，低于无雪线H_L的无云像元都为无雪；在H_S、H_L之间等间距划分高程，并且高程大于等于H_S和小于等于H_L各为一类；

然后，在每个高程带中，根据南、北和东西三类坡向继续划分高程；

在高程、坡向分类基础上，每个云像元在其所在类中的积雪概率为P_H，具体公式如下：

P_H＝N_s/N_nocloud       (1)；

其中，N_s是该类中积雪像元数目，N_nocloud是该类中非云像元数目；

P_T是当前云像元在时间上前后各5个非云像元中的积雪概率，公式如下：

$P_T=\frac{\mathrm{\Σ}1/d_{\mathrm{index}}}{\underset{i=1}{\overset{10}{\mathrm{\Σ}}}1/d_i},\mathrm{index}\∈\{\mathrm{1,2,3},...,10\}---\left(2\right);$

其中，d_i是非云像元i到当前云像元的时相距离，d_index是非云像元中积雪像元到当前云像元的时相距离。

步骤3：利用加权平均方法获得云像元的积雪概率，评判云像元是积雪还是无雪；具体实现过程是，根据下列公式得到云像元积雪的总体概率P_s，当P_s大于等于P₀时，将该云像元判断为积雪，反之为无雪：

P_s＝w₁×P_H+w₂×P_T     (3)；

其中，P₀为设定的阈值，w₁和w₂分别是P_H和P_T的权重，且w₁+w₂＝1。

所述的w₁和w₂，其计算方法是依据积雪单元中的非云像元的信息获取，具体步骤如下：

首先，在每个积雪单元中，利用步骤3中的方法求得非云像元在空间和时间上的积雪概率；

然后，在[0,1]区间内，以0.01为步长设定w₁，以公式(3)求得非云像元的积雪概率，判断其为积雪还是无雪，并将判断结果视为非云像元的预测值；

最后，将非云像元的原始值作为真实值，算出所述的遥感逐日积雪产品去云方法在非云像元上的精度oa，公式如下所示：

$\mathrm{oa}=\frac{N_s^s+N_{\mathrm{ns}}^{\mathrm{ns}}}{N_{\mathrm{nc}}}---\left(4\right);$

其中，为积雪单元中非云像元的预测值和真实值都为积雪的数目，为积雪单元中非云像元的预测值和真实值都为无雪的数目，N_nc为积雪单元中非云像元的总数目。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (5)

1.一种结合时空信息的遥感逐日积雪产品去云方法，其特征在于，包括以下步骤：

步骤1：将原始的遥感逐日积雪产品重分类为积雪、无雪和云三类，同时根据DEM数据和土地利用数据对遥感逐日积雪产品进行积雪单元划分；

步骤2：利用DEM数据和多时相数据分别计算每个云像元在空间和时间上的积雪概率；

步骤3：利用加权平均方法获得云像元的积雪概率，评判云像元是积雪还是无雪。

2.根据权利要求1所述的结合时空信息的遥感逐日积雪产品去云方法，其特征在于：步骤1中所述的对遥感逐日积雪产品进行积雪单元划分，划分后将云量大于70％的积雪单元进行合并，形成新的积雪单元。

3.根据权利要求1所述的结合时空信息的遥感逐日积雪产品去云方法，其特征在于：步骤2中所述的计算每个云像元在空间和时间上的积雪概率，具体实现过程是：

首选在每个积雪单元的高程中，找到雪线H_S和无雪线H_L；其中，高于雪线H_S的无云像元都为积雪，低于无雪线H_L的无云像元都为无雪；在H_S、H_L之间等间距划分高程，并且高程大于等于H_S和小于等于H_L各为一类；

然后，在每个高程带中，根据南、北和东西三类坡向继续划分高程；

在高程、坡向分类基础上，每个云像元在其所在类中的积雪概率为P_H，具体公式如下：

P_H＝N_s/N_nocloud             (1)；

其中，N_s是该类中积雪像元数目，N_nocloud是该类中非云像元数目；

P_T是当前云像元在时间上前后各5个非云像元中的积雪概率，公式如下：

$P_T=\frac{\mathrm{\Σ}1/d_{\mathrm{index}}}{\underset{i=1}{\overset{10}{\mathrm{\Σ}}}1/d_i},\mathrm{index}\∈\{\mathrm{1,2,3},...,10\}---\left(2\right);$

其中，d_i是非云像元i到当前云像元的时相距离，d_index是非云像元中积雪像元到当前云像元的时相距离。

4.根据权利要求3所述的结合时空信息的遥感逐日积雪产品去云方法，其特征在于：步骤3中的具体实现过程是，根据下列公式得到云像元的总体积雪概率P_s，当P_s大于等于P₀时，将该云像元判断为积雪，反之为无雪：

P_s＝w₁×P_H+w₂×P_T         (3)；

其中，P₀为设定的阈值，w₁和w₂分别是P_H和P_T的权重，且w₁+w₂＝1。

5.根据权利要求4所述的结合时空信息的遥感逐日积雪产品去云方法，其特征在于：步骤3中所述的w₁和w₂，其计算方法是依据积雪单元中的非云像元的信息获取，具体步骤如下：

首先，在每个积雪单元中，利用步骤3中的方法求得非云像元在空间和时间上的积雪概率；

然后，在[0,1]区间内，以0.01为步长设定w₁，以公式(3)求得非云像元的积雪概率，判断其为积雪还是无雪，并将判断结果视为非云像元的预测值；

最后，将非云像元的原始值作为真实值，算出所述的遥感逐日积雪产品去云方法在非云像元上的精度oa，公式如下所示：

$\mathrm{oa}=\frac{N_s^s+N_{\mathrm{ns}}^{\mathrm{ns}}}{N_{\mathrm{nc}}}---\left(4\right);$

其中，为积雪单元中非云像元的预测值和真实值都为积雪的数目，为积雪单元中非云像元的预测值和真实值都为无雪的数目，N_nc为积雪单元中非云像元的总数目。