CN115359369B - 一种基于时相自适应的山地卫星影像融合方法和系统 - Google Patents

Abstract

本发明提供了一种基于时相自适应的山地卫星影像融合方法和系统，包括，获取多个遥感数据；对多个遥感数据进行山地地形辐射校正，得到预处理遥感数据；划分得到多个时间窗口的预处理遥感数据；对于每个时间窗口中的预处理遥感数据的每一个校正像元：基于校正像元的反射率和光谱指数，判断校正像元的地表覆被类型和时间窗口内的覆被状态；基于地表覆被类型和所述覆被状态，确定筛选算法；基于筛选算法，筛选出待融合像元，并对待融合像元执行时相加权融合，得到融合像元；将每个时间窗口的融合像元进行时相加权融合，得到卫星影像数据；以消除遥感数据中云、地形等的干扰，获取时空连续的地表参数。

Description

一种基于时相自适应的山地卫星影像融合方法和系统

技术领域

本发明涉及卫星影像处理技术领域，具体而言，涉及一种基于时相自适应的山地卫星影像融合方法和系统。

背景技术

遥感影像及其估算获得的生态参量的时空连续性对于许多地表相关的应用都十分重要。近年来，国内外学者根据不同的应用目的发展了多种卫星影像时空融合方法，充分利用卫星影像的时-空-谱互补信息以克服遥感数据的时空缺失问题。然而由于云、气溶胶、传感器故障、山地地形辐射畸变等多种因素的影响，单一时相的卫星影像地表信息在空间上往往存在不连续问题，进而引起了多时相影像地表信息的时空不连续，而这种时空不连续也导致了遥感应用及遥感反演获得的地表参数时空不连续。

有鉴于此，本说明书中的一些实施例提供了一种基于时相自适应的山地卫星影像融合方法和系统，以消除遥感数据中山地地形、云等的干扰，以获取时空连续的地表参数。

发明内容

本发明的目的在于提供一种基于时相自适应的山地卫星影像融合方法，包括，按时间序列获取某一时间段的多个遥感数据；所述多个遥感数据包括同一视野；对所述多个遥感数据进行山地地形辐射校正，得到预处理遥感数据；所述预处理遥感数据为归一化至统一观测几何条件的遥感数据；对时间段进行划分，得到多个时间窗口的预处理遥感数据；对于每个时间窗口中的预处理遥感数据的每一个校正像元：基于所述校正像元的反射率和光谱指数，判断所述校正像元的地表覆被类型和时间窗口内的覆被状态；基于所述地表覆被类型和所述覆被状态，确定筛选待融合像元的筛选算法；基于所述筛选算法，从该时间窗口的预处理遥感数据中筛选出所述待融合像元，并对所述待融合像元执行时相加权融合，得到该时间窗口的融合像元；所述待融合像元为晴空像元；将每个时间窗口的融合像元进行时相加权融合，得到各时间窗口融合后的卫星影像数据。

进一步的，还包括，将所述多个遥感数据中反射率低于0.005的原像元视为噪声，将其反射率赋值为0.005。

进一步的，所述判断所述校正像元的地表覆被类型和时间窗口内的覆被状态，包括，基于所述校正像元的NDVI数据，确定所述校正像元是否为植被像元；对于所述植被像元，基于所述校正像元的NDVI数据，确定所述校正像元所在的时间窗口是否处于植被期；对于时间窗口处于非植被期的植被像元，基于所述校正像元的NDWI数据，确定所述校正像元是否处于裸地期；对于时间窗口处于非裸地期的裸地像元，基于所述校正像元的近红外波段，确定所述校正像元处于水体期还是积雪期；对于非植被像元，基于所述校正像元的NDWI数据，确定所述校正像元是否为裸地像元；对于所述裸地像元，基于所述校正像元的NDWI数据，确定所述校正像元所在的时间窗口是处于裸地期还是积雪期；对于非裸地像元，基于所述校正像元的近红外波段，确定所述校正像元是否为水体像元；对于所述水体像元，基于所述校正像元的近红外波段，确定所述校正像元所在的时间窗口是处于积雪期还是水体期。

进一步的，确定所述筛选算法，包括，对于所在的时间窗口处于植被期的所述校正像元，选择NDVI阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于植被期的所述待融合像元；对于所在的时间窗口处于裸地期的校正像元，选择NDWI阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于裸地期的所述待融合像元；对于所在的时间窗口处于积雪期或水体期的校正像元，选择短波红外波段阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于积雪期或水体期的所述待融合像元。

进一步的，所述NDVI阈值条件法为将NDVI数值大于第一阈值条件的所述校正像元作为所述待融合像元；所述NDWI阈值条件法为将NDWI数值小于第二阈值条件的所述校正像元作为所述待融合像元；所述短波红外波段阈值条件法为将短红外波段小于第三阈值条件的所述校正像元作为所述待融合像元。

进一步的，融合后的卫星影像数据包括地表反射率融合数据、融合日期数据、融合算法选择结果和各时间窗口的覆被状态数据。

进一步的，所述某一时间段的长度为一年，按时间序列获取多个遥感数据为每天获取一个遥感数据，按日期获取一年内的多个遥感数据，还包括，按天数将全年的遥感数据分为多个时段，筛选每个时段的所述待融合像元，根据所述待融合像元的时相权重进行反距离加权求和得到所述融合像元。

进一步的，当所述筛选算法无法在该时间窗口内筛选出所述待融合像元时，从相邻时间段内对应的时间窗口中筛选所述待融合像元，并标记所述待融合像元的时间。

进一步的，获取遥感数据，包括，卫星多次绕地球的经度或纬度飞行，以获取多个初始遥感数据；将所述初始遥感数据进行镶嵌拼接，得到所述遥感数据。

本发明的目的在于提供一种基于时相自适应的山地卫星影像融合系统，其特征在于，包括遥感数据获取模块、预处理模块、时间窗口划分模块、融合像元确定模块和卫星影像数据融合模块；所述遥感数据获取模块用于按时间序列获取某一时间段的多个遥感数据；所述多个遥感数据包括同一视野；所述预处理模块用于对所述多个遥感数据进行山地地形辐射校正，得到预处理遥感数据；所述预处理遥感数据为归一化至统一观测几何条件的遥感数据；所述时间窗口划分模块用于对时间段进行划分，得到多个时间窗口的预处理遥感数据；所述融合像元确定模块用于对于每个时间窗口中的预处理遥感数据的每一个校正像元：基于所述校正像元的反射率和光谱指数，判断所述校正像元的地表覆被类型和时间窗口内的覆被状态；基于所述地表覆被类型和所述覆被状态，确定筛选待融合像元的筛选算法；基于所述筛选算法，从该时间窗口的预处理遥感数据中筛选出所述待融合像元，并对所述待融合像元执行时相加权融合，得到该时间窗口的融合像元；所述待融合像元为晴空像元；所述卫星影像数据融合模块用于将每个时间窗口的融合像元进行时相加权融合，得到各时间窗口融合后的卫星影像数据。

本发明实施例的技术方案至少具有如下优点和有益效果：

本说明书中的一些实施例通过将年尺度的时间序列观测图像分为多个时间段，多个时间段以季节划分，可以在融合遥感影像的时候除去山地地表季相对影像融合的影响，提高最终影像的质量。

本说明书中的一些实施例考虑到云及其阴影、气溶胶、传感器故障、地形辐射畸变等多种因素的影响造成的时间序列卫星影像光谱辐射畸变和时空不连续问题，提出了一种基于时相自适应的山地卫星影像融合方法；该方法考虑了地表覆被状态会随着季相的变化而改变，相应的，也需要选择适用的融合规则，因此充分利用时间序列的日反射率数据的光谱和时相信息进行全局时相光谱和局部时相光谱分析，判断局部时间合成窗口内的地表覆盖状态，进而结合各合成方法的优缺点自适应的选择适用方法筛选待融合像元，并最终进行时相加权融合。该方法克服了传统直接采用云掩膜数据进行融合时云掩膜误差带来的不确定性，同时通过直接采用卫星观测分析不同地表覆盖状态的敏感波段或敏感指数时间序列特征，自适应的选择不同地表覆盖状态的最佳筛选和融合方法，避免了采用单一时相融合算法在开展大区域山地卫星影像融合时的适用性。

附图说明

图1为本发明一些实施例提供的一种基于时相自适应的山地卫星影像融合方法的示例性流程图；

图2为本发明一些实施例提供的判断像元地表覆被类型的示例性示意图；

图3为本发明一些实施例提供的一种山地卫星影像融合系统的示例性模块图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

图1为本发明一些实施例提供的一种基于时相自适应的山地卫星影像融合方法的示例性流程图。在一些实施例中，流程100可以由系统300执行。如图1所示，流程100可以包括以下步骤：

步骤110，按时间序列获取某一时间段的多个遥感数据。

在一些实施例中，某一时间段的长度可以为一年，按时间序列获取多个遥感数据为每天获取一个遥感数据，按日期获取一年内的多个遥感数据。相邻时间段可以为相邻的年份。例如，对于2018年，其相邻的时间段为2017和2019年。每天获取的遥感数据可以为全球的卫星图像。例如，极轨卫星每天至少获取全球的卫星图像。

在一些实施例中，多个遥感数据包括同一视野，同一视野可以为整个地球的观测图像，卫星可以多次绕地球的经度或纬度飞行，以获取多个初始遥感数据，将初始遥感数据进行拼接，得到遥感数据。

例如，极轨卫星围绕地球从北极向南极飞行，在每条轨道上获取地球的一部分观测图像，将地球多个部分的观测图像进行拼接，可以得到全球的地表反射率数据。

步骤120，对多个遥感数据进行山地地形辐射校正，得到预处理遥感数据。

预处理遥感数据为归一化至统一观测几何条件的遥感数据。例如，采用SCS+C地形辐射校正模型将多个山地遥感数据校正为水平地表条件下的拥有统一观测几何条件的地表反射率。

在一些实施例中，还可以将多个遥感数据中反射率低于0.005的原像元视为噪声，将其反射率赋值为0.005。原像元为获取的遥感数据的像元。对于经过预处理后的遥感数据，其像元可以被称为校正像元。

步骤130，对时间段进行划分，得到多个时间窗口的预处理遥感数据。

在一些实施例中，遥感数据可以为每天的遥感数据，由对遥感数据进行预处理得到的预处理遥感数据可以为每天的预处理遥感数据，时间段可以为一年，每个时间窗口可以包括多天的预处理遥感数据。例如，每个时间窗口可以包括16天的预处理遥感数据，对于一年的预处理遥感数据可以有23个时间窗口。当然，每个时间窗口还可以包括其他天数的预处理遥感数据，可以根据实际需求设置。

多个预处理遥感数据包括同一视野，对于视野中的位置，每个预处理遥感数据中有与之相对的校正像元，其中，处理每个时间窗口中的预处理遥感数据的每一个校正像元可以是指处理同一视野在多个预处理遥感数据中的每一个校正像元。例如，对于视野中的位置A，和预处理遥感数据1、预处理遥感数据2和预处理遥感数据3，其中，校正像元A1是预处理遥感数据1中获取的位置A的像元，同理A2和A3，可以对A1、A2和A3进行融合处理，得到能代表位置A处的实际观测数据的融合像元，该像元可以被认为是最佳像元。

得到融合像元的方法可以包括以下内容：

步骤140，基于校正像元的反射率和光谱指数，判断校正像元的地表覆被类型和时间窗口内的覆被状态。

地表覆被类型可以用于反映地表的情况。地表覆被类型包括永久冰川积雪、裸地（裸土、裸岩）、永久水体、水体（冬季冻结）、热带常绿植被、温带落叶植被（落叶林、农田、草地等）、寒带落叶植被（冬季被雪覆盖）等七种。

不同地表覆被类型对不同波段的反射率和光谱指数不同，因此，可以基于地表对敏感波段的反射率和光谱指数确定地表覆被类型。针对以上七种地表覆被类型，可以选择可见光蓝波段、红波段、近红外波段和短波红外波段等波段作为敏感波段，并基于敏感波段的反射率构建归一化植被指数（NDVI）、归一化水体指数（NDWI）和归一化积雪指数（NDSI）作为判断地表状态的敏感光谱。

时间窗口内的覆被状态可以是指在该时间段内，地表的覆被类型。例如，像元B，其在夏季的地表覆被类型可以为草地，其在冬季的地表覆被类型可以为积雪。

在一些实施例中，可以基于不同地表覆被类型对敏感波段的反射率和光谱指数，确定地表覆被类型。例如，通过将待测的像元的反射率和光谱指数与典型值进行比对，即可得到待测像元处的地表覆被类型。

典型值可以为分析典型数据得到的不同地表覆被类型在有云和无云状态下的标准反射率和标准光谱指数。例如，在VIIRS全球地表反射率数据中，在全球土地覆被数据的辅助下，利用目视解译的方法从1月和6月中两个云覆盖较少的日期的地表反射率影像中选择了植被、裸地、冰雪、水体、厚云的纯像元各约2×10⁵个，分析其敏感波段的光谱分布规律，得到在无云干扰的条件下的典型值。同理，在VIIRS全球日地表反射率数据中选择永久冰川积雪、裸地、永久水体、水体（冬季冻结）、热带常绿植被、温带落叶植被（落叶林、农田、草地等）、寒带落叶植被（冬季被雪覆盖）七种不同地表类型的全年时间序列数据，分析其各类型的全年反射率数据和光谱指数在有云干扰的条件下的分布规律，得到在有云干扰的条件下的典型值。其中，为了避免因冬季降雪或冻结的影像，永久水体像元主要从热带区域（10°N-10°S）选择，冻结水体主要在北半球寒带区域（纬度>40°N）选择，寒带落叶植被主要在北半球寒带区域（纬度>40°N）选择，以包含其冬季积雪像元。

在一些实施例中，可以基于年尺度，确定像元全年内的地表覆被类型。

例如，通过分析一年内所有时间序列的像元光谱频率特征进行地表覆盖状态的判定。首先根据像元全年时间序列的光谱反射率或光谱指数累积特征进行分析，年尺度地表覆被状态YCT具体判别为：

（1）

式中，N_observation是满足括号内条件假设的全年有效观测数量，R_yearly表示选择的敏感光谱波段或光谱指数，T_R表示敏感光谱波段或光谱指数的阈值，T_N表示有效观测数量的阈值。该判别规则的目的是判别地表至少在一定时间内是否有植被覆盖，否则年内地表为非植被，主要被冰雪、水体或裸露地表覆盖。

然后，根据地表覆盖数据，利用分层抽样法选择不同土地覆盖类型的样本，并分为两类：有显著季节变化和无显著季节变化样本。其中无显著季节变化覆被类型样本包括永久水体、永久积雪、裸地（裸土、裸岩）和常绿植被。有显著季节变化类型包括冬季结冰水体、落叶植被（如落叶林、灌木、草本植被和农田）。对于每一类覆被类型样本，选择包含有云和无云的时间序列卫星观测样本，分别分析所选敏感波段和敏感光谱指数的直方图和累积观测频率，选择累积观测频率大于95%处的光谱或指数阈值作为阈值条件。

在另一些实施例中，可以基于局部时间窗口的变化特征，确定像元时段内的地表覆被类型。局部时间窗口可以为一段时间，例如，16天。

例如，对于永久水体或永久积雪覆盖，本发明选择单一最佳适宜方法进行待融合像元筛选。对地表覆被状态会随季相变化的非永久水体和永久积雪类型，其地表覆被类型状态TCT根据局地时间窗口内观测的累积观测频率进行确定：

（2）

其中，N_{T-observation}表示局部时间窗口内的累积有效观测数量，T_TN代表该时间窗口内满足条件的观测数量的阈值。R_sub-yearly表示选择的敏感光谱波段或光谱指数，该判别规则的目的是判别地表至少在该时间窗口的一定时间内是否为水体，否则时间窗口内地表积雪。

关于确定校正像元的地表覆被类型的更多内容，参见图2及其相关描述。

步骤150，基于地表覆被类型和覆被状态，确定筛选待融合像元的筛选算法。

在一些实施例中，可以基于校正像元的光谱指数对校正像元进行筛选。在另一些实施例中，可以将相近的地表对应的校正像元作为一组像元，然后对每组像元进行融合。关于确定对校正像元进行融合的融合算法的更多内容，参见图2及其相关描述。

步骤160，基于筛选算法，从该时间窗口的预处理遥感数据中筛选出待融合像元，并对待融合像元执行时相加权融合，得到该时间窗口的融合像元；待融合像元为晴空像元。

多个待融合像元可以为同一个时间窗口无云的观测时间。例如，可以将多个时序的遥感数据中均能反映裸地状态的遥感影像对应的获取时间作为观测裸地的待融合时相。融合像元可以是将多个待融合时相进行加权求和，以便反映该位置的地表状况。关于待融合像元的更多内容，参见图2及其相关描述。

在一些实施例中，在获得同一时间窗口多个待融合时相观测数据后，融合像元可以根据各待融合像元日序与目标融合像元日序的差值进行日序反距离加权求和确定，其中，目标融合像元日序可以为时间窗口的中值：

（3）

其中，R _fusion 表示融合像元的反射率，R _i 为第i个待融合像元的反射率，N为待融合像元的个数，w _i为第i个待融合像元权重，根据第i个待融合像元日序DOY _i 与目标融合像元日序DOY _fusion 差值进行确定：

（4）

在一些实施例中，由于地表反射率数据中可能会存在缺失数据和噪声数据，因此，当各波段反射率低于0.005且非填充值时，由于传感器的探测能力一般无法探测此类弱反射率信息，因此，将该波段反射率赋值为0.005，以对数据进行降噪处理。此外，由于山地局地地形起伏，导致卫星和太阳以及目标地物的观测几何关系发生改变，进而导致同一地表覆被类型其因观测几何关系的改变而地表反射率发生变化。因此本发明进一步利用STS+C地形辐射校正模型，将山地地表反射率归一化至平坦地表统一观测几何条件，增加用于融合的多时相遥感数据的可比性和一致性。

在一些实施例中，当筛选算法无法在该时间窗口内筛选出待融合像元时，从相邻时间段内对应的时间窗口中筛选待融合像元，并标记待融合像元的时间。

例如，时间窗口均无满足筛选算法条件的云像元时，采用相邻年份对应时间窗口数据进行融合并标记时间，包括，

对于时间窗口内无满足阈值条件的待融合像元，检查待融合时间窗口相邻年份的对应时间窗口内满足阈值条件的像元，采用时相加权方法进行融合，并标记融合像元的时段。

在一些实施例中，一年的部分时间窗口内的观测可能全部受云干扰，在这种情况下，将无法获得满足筛选算法的待融合像元。此时对于视野中的区域，检索该区域相邻前后1年的遥感数据，针对相同的时间窗口判断是否存在满足筛选条件的待融合数据，按照时相加权算法进行融合，并标记该融合像元来自于前后1年的数据时间。

步骤170，将每个时间窗口的融合像元进行时相加权融合，得到各时间窗口融合后的卫星影像数据。

在一些实施例中，可以将获取的视野中所有位置对应的待融合像元构成的影像数据进行时相加权，作为该多个遥感数据融合后的影像数据。其中，对于视野中的多个像元，可能来自不同时相。例如，可以采用逐行读写的方式进行数据行读取，获取每个像元年逐日的时间序列反射率和光谱指数，然后对数据进行降噪处理和地形辐射校正后进行全年和局部时间窗口内地表覆被状态的判断和筛选算法的自适应选择，基于选择的筛选算法获取待融合的遥感数据，最后输出数据包括融合结果数据层、融合日期数据层、融合算法数据层、以及融合窗口内地表覆被状态数据层。数据采用ENVI标准格式，采用正弦投影方式进行存储。

在一些实施例中，融合后的卫星影像数据可以包括12个数据层，分别是9个波段的地表反射率合成数据，1个波段的合成日期数据、1个波段的筛选算法数据和1个波段的16天地表状态数据。数据采用ENVI标准格式，采用正弦投影方式进行存储。

在一些实施例中，遥感数据可以为卫星一年或多年的监测数据，可以按天数将全年的遥感数据分为多个时段（8天、10天、15天、月等），然后，计算每个时段中的待融合像元，并通过时相加权获取该时段的融合像元；基于每个时段的融合像元，得到多个时段的融合后的卫星影像数据。

图2为本发明一些实施例提供的判断像元地表覆被类型的示例性示意图。在一些实施例中，对于所在的时间窗口处于植被期的校正像元，选择NDVI阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于植被期的待融合像元；对于所在的时间窗口处于裸地期的校正像元，选择NDWI阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于裸地期的待融合像元；对于所在的时间窗口处于积雪期或水体期的校正像元，选择短波红外波段阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于积雪期或水体期的待融合像元。在一些实施例中，流程200可以由融合像元确定模块340执行。如图2所示，流程200可以包括以下步骤：

步骤210，基于校正像元的NDVI数据，确定校正像元是否为植被像元。

例如，在当像元全年NDVI数据小于0.2的天数大于全年有效观测天数的95%时，算法将该像元判别为非植被像元，否则为植被像元。

步骤220，对于植被像元，基于校正像元的NDVI数据，确定校正像元所在的时间窗口是否处于植被期。

例如，对于植被像元，直接进行局部时间窗口的判断，若时间窗口内的NDVI小于0.2的天数大于有效观测的95%，则在该时间窗口内该像元主要为非植被期（或全部为云覆盖）。

步骤230，对于时间窗口处于非植被期的植被像元，基于校正像元的NDWI数据，确定校正像元是否处于裸地期。

例如，对于非植被期的像元，算法进一步利用NDWI进行判断，若时间窗口内NDWI大于-0.115的天数大于90%，则表明该时间窗口内没有裸地期的像元，否则该时间窗口内该像元被裸地和云覆盖。

步骤240，对于时间窗口处于非裸地期的裸地像元，基于校正像元的近红外波段，确定校正像元处于水体期还是积雪期。

例如，对于非裸地期的像元，其有可能为水体（季节性淹水湿地）和积雪两种情况，进一步利用近红外波段中积雪和水体的光谱差异进行判断，若时间窗口内近红外波段大于0.2的反射率天数大于有效观测天数90%，则在该时间窗口内该像元主要为积雪和云（或全部为云），否则为水体。

步骤250，对于非植被像元，基于校正像元的NDWI数据，确定校正像元是否为裸地像元。

例如，对于非植被像元，由于裸地与水体、积雪相比拥有完全不同的NDWI指数特征，即裸地的NDWI为负值而其他非植被地表类型为正值，因此，选择NDWI进行下一层次的判断。当非植被像元全年NDWI大于-0.115（水体的NDWI下限）的天数大于全年有效观测的95%时，则认为不包含裸地像元；否则该像元为裸地（可能季节性被积雪覆盖）。

步骤260，对于裸地像元，基于校正像元的NDWI数据，确定校正像元所在的时间窗口是处于裸地期还是积雪期。

例如，对于裸地像元，进一步利用NDWI在时间窗口内对地表进行判断，若时间局部窗口内的NDWI大于-0.115的像元超过该时间窗口内有效观测的95%，则该时间窗口内为积雪和云（或者全部为云），否则为裸地和云。

步骤270，对于非裸地像元，基于校正像元的近红外波段，确定校正像元是否为水体像元。

例如，对于非裸地像元，主要进行积雪和水体的判别，选择在二者光谱反射率有显著差异近红外波段作为区分。若非裸地像元中近红外波段地表反射率大于0.2的全年有超过有效观测的90%，则认为该时间窗口内主要为永久积雪；否则为水体。

步骤280，对于水体像元，基于校正像元的近红外波段，确定校正像元所在的时间窗口是处于积雪期还是水体期。

由于水体像元可能受冬季冰雪覆盖。进一步利用近红外波段反射率进行水体和冰雪的区分。

在一些实施例中，NDVI阈值条件法为将NDVI数值大于第一阈值条件的校正像元作为待融合像元；NDWI阈值条件法为将NDWI数值小于第二阈值条件的校正像元作为待融合像元；短波红外波段阈值条件法为将短红外波段小于第三阈值条件的校正像元作为待融合像元。其中，第一阈值条件、第二阈值条件和第三阈值条件为预先设置的筛选条件。通过将视野中的待融合像元进行时相加权输出，得到观测清晰的融合后的遥感数据。例如，确定对校正像元进行融合的筛选算法，包括，对于处于植被期的像元，选择NDVI值大于0.2的像元作为待融合像元；对于处于裸地期的像元，选择NDWI值小于-0.115的像元作为待融合像元；对于处于积雪期和水体期的像元，选择短波红外波段反射率小于0.2的像元作为待融合像元。

例如，在地表覆盖状态判定后，基于地表覆盖类型选择筛选算法进行待融合像元的选择。其中，当局部时间窗口内有植被或裸土观测时选择归一化植被指数大于阈值条件进行待融合像元选择；当局部时间窗口内主要为水体或冰雪时，由于短波红外波段水体、冰雪的反射率均比云低，可选择短波红外反射率小于阈值条件进行该时间窗口内无云像元的选择。最佳筛选算法自适应选择的算法为：

（5）

其中，

代表时间窗口T内经算法筛选获得的无云观测像元（待融合 像元），NDVI代表归一化植被指数，T代表融合时间窗口大小，单位为天，

代表短波 红外波段反射率，TCT=Veg or Barren 代表合成时间窗口T内的地表覆盖状态为植被或裸 地，TCT=Water or Snow/ice代表合成时间窗口T内的地表覆盖状态为水体或冰雪，n为待融 合像元（无云观测像元）的个数。

图3为本发明一些实施例提供的一种山地卫星影像融合系统的示例性模块图。如图3所示，系统300可以包括遥感数据获取模块310、预处理模块320、时间窗口划分模块330、融合像元确定模块340和卫星影像数据融合模块350。

遥感数据获取模块310用于按时间序列获取某一时间段的多个遥感数据；多个遥感数据包括同一视野。关于遥感数据获取模块310的更多内容，参见图1及其相关描述。

预处理模块320用于对多个遥感数据进行山地地形辐射校正，得到预处理遥感数据；预处理遥感数据为归一化至统一观测几何条件的遥感数据。关于预处理模块320的更多内容，参见图1及其相关描述。

时间窗口划分模块330用于对时间段进行划分，得到多个时间窗口的预处理遥感数据。关于时间窗口划分模块330的更多内容，参见图1及其相关描述。

融合像元确定模块340用于对于每个时间窗口中的预处理遥感数据的每一个校正像元：基于校正像元的反射率和光谱指数，判断校正像元的地表覆被类型和时间窗口内的覆被状态；基于地表覆被类型和覆被状态，确定筛选待融合像元的筛选算法；基于筛选算法，从该时间窗口的预处理遥感数据中筛选出待融合像元，并对待融合像元执行时相加权融合，得到该时间窗口的融合像元；待融合像元为晴空像元。关于融合像元确定模块340的更多内容，参见图1及其相关描述。

卫星影像数据融合模块350用于将每个时间窗口的融合像元进行时相加权融合，得到各时间窗口融合后的卫星影像数据。关于卫星影像数据融合模块350的更多内容，参见图1及其相关描述。

通过将NASA NPP 8天合成的地表反射率产品、MODIS 标准的8天地表反射率合成产品MYD09A1 （C6）与通过本说明书中一些实施例的方法获得的8天地表反射率数据进行对比：

在云筛选方面，本发明算法尽管未采用云掩膜数据，其融合影像在残云量上要明显优于其余两种产品。NASA NPP 8天合成产品与MODIS的MYD09标准合成产品均采用了云掩膜产品作为辅助数据。在合成影像中，两种产品均出现很多噪声点及不规则边缘，且条带较为明显。而本发明算法获得的合成结果在影像总体质量、纹理特征和清晰度上要优于MODIS影像合成产品。其中对比区域为从10°N至60°N的纬度梯度，基本涵盖了从热带到寒带的气候带，且覆盖青藏高原，具有很强的代表性。对比的合成日期为105-112天。

对于生成的图像，本发明合成方法的时间序列曲线具有更强的稳定性和光谱一致性。具体而言，对永久冰川积雪，本发明选择的合成结果其短波红外波段较低且整个时间序列上一致性较好，这与积雪在该波段的吸收特性较强是一致的。其中，NASA VIIRS的合成产品反射率波动较为剧烈。对水体而言，由于水体在近红外波段具有更强的吸收特性，因此其NDVI往往为负值。而NASA VIIRS的合成产品部分合成日期的NDVI接近于零，表明其合成影像未能正确选择水体像元。对植被而言，本发明合成方法能够更好的选择植被的上包络线，表示该合成期内植被的最佳生长状态。而对于沙漠而言， NASA VIIRS的合成影像中仍存在少量残云或云阴影的干扰。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于时相自适应的山地卫星影像融合方法，其特征在于，包括，

按时间序列获取某一时间段的多个遥感数据；所述多个遥感数据包括同一视野；

对所述多个遥感数据进行山地地形辐射校正，得到预处理遥感数据；所述预处理遥感数据为归一化至统一观测几何条件的遥感数据；

对时间段进行划分，得到多个时间窗口的预处理遥感数据；

对于每个时间窗口中的预处理遥感数据的每一个校正像元：

基于所述校正像元的反射率和光谱指数，判断所述校正像元的地表覆被类型和时间窗口内的覆被状态；

基于所述地表覆被类型和所述覆被状态，确定筛选待融合像元的筛选算法；所述筛选算法为基于所述校正像元的光谱指数对所述校正像元进行筛选，所述筛选算法的表达式为：

其中，

代表时间窗口T内经算法筛选获得的无云观测像元，即，待融合像 元，NDVI代表归一化植被指数，T代表融合时间窗口大小，单位为天，

代表短波红 外波段反射率，TCT=Veg or Barren 代表合成时间窗口T内的地表覆盖状态为植被或裸地， TCT=Water or Snow/ice代表合成时间窗口T内的地表覆盖状态为水体或冰雪，n为待融合 像元的个数；

基于所述筛选算法，从该时间窗口的预处理遥感数据中筛选出所述待融合像元，并对所述待融合像元执行时相加权融合，得到该时间窗口的融合像元；所述待融合像元为晴空像元；

将每个时间窗口的融合像元进行时相加权融合，得到各时间窗口融合后的卫星影像数据。

2.根据权利要求1所述的基于时相自适应的山地卫星影像融合方法，其特征在于，还包括，将所述多个遥感数据中反射率低于0.005的原像元视为噪声，将其反射率赋值为0.005。

3.根据权利要求1所述的基于时相自适应的山地卫星影像融合方法，其特征在于，所述判断所述校正像元的地表覆被类型和时间窗口内的覆被状态，包括，

基于所述校正像元的NDVI数据，确定所述校正像元是否为植被像元；

对于所述植被像元，基于所述校正像元的NDVI数据，确定所述校正像元所在的时间窗口是否处于植被期；

对于时间窗口处于非植被期的植被像元，基于所述校正像元的NDWI数据，确定所述校正像元是否处于裸地期；

对于时间窗口处于非裸地期的裸地像元，基于所述校正像元的近红外波段，确定所述校正像元处于水体期还是积雪期；

对于非植被像元，基于所述校正像元的NDWI数据，确定所述校正像元是否为裸地像元；

对于所述裸地像元，基于所述校正像元的NDWI数据，确定所述校正像元所在的时间窗口是处于裸地期还是积雪期；

对于非裸地像元，基于所述校正像元的近红外波段，确定所述校正像元是否为水体像元；

对于所述水体像元，基于所述校正像元的近红外波段，确定所述校正像元所在的时间窗口是处于积雪期还是水体期。

4.根据权利要求3所述的基于时相自适应的山地卫星影像融合方法，其特征在于，确定所述筛选算法，包括，

对于所在的时间窗口处于植被期的所述校正像元，选择NDVI阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于植被期的所述待融合像元；

对于所在的时间窗口处于裸地期的校正像元，选择NDWI阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于裸地期的所述待融合像元；

对于所在的时间窗口处于积雪期或水体期的校正像元，选择短波红外波段阈值条件法从该时间窗口对应的预处理遥感数据中筛选出多个处于积雪期或水体期的所述待融合像元。

5.根据权利要求4所述的基于时相自适应的山地卫星影像融合方法，其特征在于，所述NDVI阈值条件法为将NDVI数值大于第一阈值条件的所述校正像元作为所述待融合像元；所述NDWI阈值条件法为将NDWI数值小于第二阈值条件的所述校正像元作为所述待融合像元；所述短波红外波段阈值条件法为将短红外波段小于第三阈值条件的所述校正像元作为所述待融合像元。

6.根据权利要求1所述的基于时相自适应的山地卫星影像融合方法，其特征在于，融合后的卫星影像数据包括地表反射率融合数据、融合日期数据、融合算法选择结果和各时间窗口的覆被状态数据。

7.根据权利要求1所述的基于时相自适应的山地卫星影像融合方法，其特征在于，所述某一时间段的长度为一年，按时间序列获取多个遥感数据为每天获取一个遥感数据，按日期获取一年内的多个遥感数据，还包括，

按天数将全年的遥感数据分为多个时段，筛选每个时段的所述待融合像元，根据所述待融合像元的时相权重进行反距离加权求和得到所述融合像元。

8.根据权利要求1所述的基于时相自适应的山地卫星影像融合方法，其特征在于，当所述筛选算法无法在该时间窗口内筛选出所述待融合像元时，从相邻时间段内对应的时间窗口中筛选所述待融合像元，并标记所述待融合像元的时间。

9.根据权利要求1所述的基于时相自适应的山地卫星影像融合方法，其特征在于，获取遥感数据，包括，

卫星多次绕地球的经度或纬度飞行，以获取多个初始遥感数据；

将所述初始遥感数据进行镶嵌拼接，得到所述遥感数据。

10.一种基于时相自适应的山地卫星影像融合系统，其特征在于，包括遥感数据获取模块、预处理模块、时间窗口划分模块、融合像元确定模块和卫星影像数据融合模块；

所述遥感数据获取模块用于按时间序列获取某一时间段的多个遥感数据；所述多个遥感数据包括同一视野；

所述预处理模块用于对所述多个遥感数据进行山地地形辐射校正，得到预处理遥感数据；所述预处理遥感数据为归一化至统一观测几何条件的遥感数据；

所述时间窗口划分模块用于对时间段进行划分，得到多个时间窗口的预处理遥感数据；

所述融合像元确定模块用于对于每个时间窗口中的预处理遥感数据的每一个校正像元：基于所述校正像元的反射率和光谱指数，判断所述校正像元的地表覆被类型和时间窗口内的覆被状态；基于所述地表覆被类型和所述覆被状态，确定筛选待融合像元的筛选算法；基于所述筛选算法，从该时间窗口的预处理遥感数据中筛选出所述待融合像元，并对所述待融合像元执行时相加权融合，得到该时间窗口的融合像元；所述待融合像元为晴空像元；所述筛选算法为基于所述校正像元的光谱指数对所述校正像元进行筛选，所述筛选算法的表达式为：

其中，

代表时间窗口T内经算法筛选获得的无云观测像元，即，待融合像 元，NDVI代表归一化植被指数，T代表融合时间窗口大小，单位为天，

代表短波红 外波段反射率，TCT=Veg or Barren 代表合成时间窗口T内的地表覆盖状态为植被或裸地， TCT=Water or Snow/ice代表合成时间窗口T内的地表覆盖状态为水体或冰雪，n为待融合 像元的个数；

所述卫星影像数据融合模块用于将每个时间窗口的融合像元进行时相加权融合，得到各时间窗口融合后的卫星影像数据。

Images