CN105046648B - 一种构建高时空遥感数据的方法 - Google Patents

Abstract

本发明公开了一种构建高时空遥感数据的方法，结合像元分解降尺度方法和STARFM模型的CDSTARFM方法，首先利用像元分解降尺度算法对低分辨率数据进行降尺度处理，然后用降尺度数据取代STARFM算法中直接重采样的低分辨率数据，最后利用两者相结合的CDSTARFM算法进行数据融合。该方法可以有效解决像元分解降尺度融合数据的“图斑”现象和STARFM模型寻找纯净MODIS像元难的问题。利用本发明提供的方法融合Landsat8和MODIS数据构建高时空分辨率遥感数据。结果表明，本发明提供的方法比STARFM和像元分解降尺度算法具有更高的融合精度。

Description

一种构建高时空遥感数据的方法

技术领域

本发明涉及一种构建高时空遥感数据的方法，尤其涉及一种在地块破碎区域构建高时空遥感数据的方法，属于遥感影像数据处理领域。

背景技术

高空间、时间分辨率的遥感数据在监测土地覆盖变化和作物生长以及识别地物类型等方面具有重要的作用。然而，目前获取到的卫星传感器数据存在空间分辨率、时间分辨率、光谱分辨率相互矛盾，即同时满足高空间、高时间和高光谱分辨率的卫星遥感数据是不现实的。比如，Landsat卫星的多光谱影像空间分辨率是30米，其适中的空间分辨率和易获取性在植被指数提取、监测土地覆盖动态变化和生态系统方面具有广泛的应用，但因其16天的重访周期和阴雨天气的影响，难以保证获得连续有效的影像进行农作物生长、物候变化和自然灾害的监测；搭载在Terra/Aqua卫星上的中分辨率成像光谱仪(MODIS)能提供时间分辨率为1天、空间分辨率为250-1000米的遥感影像，被广泛应用于大尺度区域的地表覆盖类型监测，但MODIS数据的低空间分辨率限制其在景观破碎和异质性较强区域的应用。若能综合中高分辨率影像(如Landsat OLI)的高分辨率特性和MODIS时间分辨率的优势，将会大幅度提高植被遥感监测能力。

近些年，国内外学者提出了系列获取高时空分辨率数据的遥感数据融合模型。这些融合模型总体分为两大类，一类是像元分解降尺度算法(Downscaling mixed pixel)，该算法是基于线性光谱混合模型，通过从高分辨率数据中提取端 元和丰度来分解低分辨率的像元获得端元光谱值。该算法简单易行，不需要中、高分数据具有相同的波段，只需要高分辨率的地表覆盖类型数据就可以进行数据融合，但是该算法获得的端元反射率是整个区域或局部区域的平均反射率，不能很好地体现地物光谱的空间差异性会出现“图斑”的现象，另一类是以Gao等提出的自适应遥感图像时空融合方法为代表，该模型用于融合Landsat和MODIS数据构建每日的Landsat数据，在融合过程中不仅考虑了空间的差异性，而且也考虑了时间的差异性，是目前应用最广泛的时空融合模型之一。但STARFM模型存在两个问题：其一，若短暂或突变的地表变化信息没有被基期的Landsat影像记录下，那么融合的影像就不能捕捉到该地表变化信息；其二，该模型构建的数据质量依赖于从MODIS数据中提取的纯净像元的时间变化信息，若在窗口内没有纯净的MODIS像元，数据融合质量在一定程度上会降低。为了解决第一个不足，Hilker等基于STARFM模型提出了一种新的融合算法。STAARCH算法分别从Landsat和MODIS数据中提取空间的变化和时间变化，通过选取最佳的基准期Landsat影像来提高融合算法的精度。对于第二问题，由于MODIS影像空间分辨率低，很难在地块破碎、异质性强的区域内找到纯净的MODIS像元，为了解决纯净MODIS像元对STARFM算法精度的影响，Zhu等提出了适合异质性区域的增强型STARFM(ESTARFM)模型，ESTARFM模型通过假设一段时间内地物反射率为线性变化以及混合像元由地物端元线性组合，在模型中引入一个转换系数来提高破碎区域数据融合的精度，但ESTARFMM模型至少需要两对同期的Landsat和MODIS基期影像，与STARFM相比增加了数据获取的难度。

发明内容

本发明针对STARFM模型因在破碎区域难搜寻到纯净像元而降低融合数据精度的问题，首先利用像元分解降尺度算法对低分辨率数据进行降尺度处 理，然后将降尺度数据取代STARFM算法中直接重采样的低分辨率数据，最后再利用两者相结合的CDSTARFM算法进行数据融合。

具体地，本发明提出了一种结合像元分解降尺度算法和STARFM模型的CDSTARFM(Combining Downscale Mixed Pixel with Spatial and Temporal AdaptiveReflectance Fusion Model)算法并用于融合Landsat8和MODIS反射率数据构建每天的Landsat8反射率数据实验。

本发明提供一种构建高时空遥感数据的方法，主要包括以下步骤：

步骤一、选取一个基期(t_k)的中、低分辨率遥感影像和时序低分辨率影像；

步骤二、对t_k时期中分辨率影像进行聚类，得到分类图像后进行类别丰度提取，得到类别丰度图；

步骤三、通过像元分解模型利用步骤二得到的类别丰度图以及t₀和t_k两期低分辨率影像进行像元分解，得到t₀和t_k两期的降尺度影像；

步骤四、利用t_k时期中分辨率影像筛选相似像元，得到窗口内中心像元的相似性像元；

步骤五、利用步骤三得到的t₀和t_k两期低分辨率影像降尺度的影像数据和步骤四得到的相似性像元，通过加权窗口内相似性像元计算预测时期的中心像元值，从而得到预测期t₀的中分辨率遥感影像。

优选的，上述步骤一中的一个基期(t_k)包括该时期的一景中分辨率遥感影像和一景低分辨率遥感影像。

优选的，上述步骤二中的丰度提取具体为通过公式(1)提取每个混合像元内各类别的丰度：

式(1)中，f_c(i，c)表示i位置低分辨率像元内c类端元的丰度，Q表示低分辨率像元内c端元的像元数，S表示低分辨率像元内所有端元的像元数。

优选的，上述步骤三具体包括：

步骤3.1混合像元的光谱值等于该像元内各地类的光谱值与其丰度的线性组合，可利用公式(2)表示。

步骤3.2选择合适窗口进行混合像元分解，利用公式(3)求解窗口内各类别的光谱值。

公式的约束条件：

其中，式(2)和式(3)中，R(i，t)为t时期i位置低分辨率像元的反射率；f_c(i，c)为i位置低分辨率像元内类别c占该像元的面积比；为t时期类别c的平均反射率；ξ(i，t)为残差；k为窗口内类别数；

步骤3.3把求得的各类别光谱值依照类别对应到窗口内相应地类的像元上，获得降尺度的数据。

优选的，上述步骤四和步骤五通过在一定窗口内筛选与中心像元光谱相似的相似性像元，利用相似性像元的光谱差异性，时间差异性和与中心像元的相对距离来加权计算预测期的中心像元值。

优选的，上述步骤四和步骤五具体包括：

步骤4.1由t_k时期中分辨率影像的green、red和NIR波段的阈值θ_b确定窗口内中心像元的相似性像元，阈值θ_b通过公式(4)计算，窗口内与中心像元差值的绝对值小于各波段θ_b的像元作为相似性像元。

式(4)中，θ_b表示窗口内b波段的阈值，N为窗口内中分辨率像元个数， x_i为i位置的像元反射率，μ为窗口内像元反射率均值，m为类别数；

步骤4.2计算相似像元的权重大小W_ijk，权重大小W_ijk有三个指标来衡量：光谱的差异性S_ijk，时间的差异性T_ijk和相似性像元与中心像元的相对距离D_ijk；

T_ijk＝|M(x_i，y_j，t_k)-M(x_i，y_j，t₀)| (6)

步骤4.3利用公式(9)加权窗口内相似像元的反射率计算预测期的中心像元反射率，生成最终的融合影像

上述公式中，w为窗口大小，L(x_i，y_j，t_k)和M(x_i，y_j，t_k)分别为t_k时期给定位置(x_i，y_j)的中分辨率数据和降尺度数据的像元值，(x_w/2，y_w/2)为窗口的中心像元，光谱差异性S_ijk值越小说明给定位置与邻近像元的光谱相似度越高，赋予较高的权重；时间差异性T_ijk值越小说明该段时间内光谱变化越小，赋予较高的权重；相对距离D_ijk值越小，赋予较高的权重。

本发明提供的构建高时空遥感数据的方法可以有效解决像元分解降尺度融合数据的“图斑”现象和STARFM模型寻找纯净MODIS像元难的问题。利用本发明提供的方法融合Landsat8和MODIS数据构建高时空分辨率遥感数据。结果表明，本发明提供的方法比STARFM和像元分解降尺度算法具有更高的融合精度及更接近真实影像的目视效果，如NIR波段，相关系数r(0.96vs0.95vs0.90；RMSE(0.0245vs 0.0300vs 0.0401)；ERGAS(0.5416vs0.6507vs 0.8737)。

附图说明

图1为本发明流程示意图；

图2为三种不同方法预测的反射率与真实反射率的散点图；

图3(a)-3(j)表示不同方法在最佳窗口下的融合影像与真实影像比较(NIR-Red-Green组合)图；其中图3(a)、图3(b)分别表示2014-08-19和09-04的Landsat8影像；图3(c)-图3(e)分别表示降尺度方法、STARFM和CDSTARFM在最佳窗口下的融合影像；图3(f)-图3(j)分别表示对应的子区域影像。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实施方式对本发明作进一步的详细描述。

本发明结合像元分解降尺度算法和STARFM模型的CDSTARFM(CombiningDownscale Mixed Pixel with Spatial and Temporal Adaptive Reflectance FusionModel)方法并用于融合Landsat8和MODIS反射率数据构建每天的Landsat8反射率数据实验。

本发明结合了像元分解降尺度算法和时空自适应融合模型STARFM的各自优势，既可以解决像元分解降尺度融合数据的“图斑”现象又可以在一定程度上增加了STARFM算法搜寻到纯净像元的概率，从而提高数据融合的精度。CDSTARFM方法输入的数据包括一个基期(t_k)的中、低分辨率遥感影像和时序低分辨率影像，输出结果是具有中分辨率特征的时序影像。CDSTARFM方法流程见图1。

一、类别定义及丰度提取

混合像元内的地物类别及其丰度是进行混合像元分解的关键。类别是通过传统方法聚类中分辨遥感影像获取，混合像元内各类别的丰度是在忽略不同传 感器数据之间的地理配准误差的前提下通过公式(1)提取的，并且类别和丰度假设在一定时间内不随时间发生变化。利用传统非监督分类方法ISO-DATA对中分辨率遥感影像聚类获得分类图像，利用ArcGIS在分类图上创建大小为低分辨率像元大小尺度的网格提取格网内各类别的丰度。

式(1)中，f_c(i，c)表示i位置低分辨率像元内c类端元的丰度，Q表示低分辨率像元内c端元的像元数，S表示低分辨率像元内所有端元的像元数。

二、降尺度数据获取

混合像元分解的基础是线性光谱混合模型，在忽略不同传感器数据之间的地理配准误差和大气校正误差的条件下，混合像元的光谱值等于该像元内端元的光谱值与其丰度的线性组合(公式2)。然而，由于地物存在空间差异性，公式(2)在一定大小的矩形窗口内进行分解混合像元，这种方式既可以消除中、高分辨率影像之间因配准误差带来的影响，又可以体现同类地物的空间差异性，但过大的窗口会削弱空间的异质性，也不一定能提高混合像元分解的精度。因此，选择合适的矩形窗口进行混合像元分解是必要的。为选择一个合适窗口，设置不同大小的窗口进行混合像元分解，对各窗口下降尺度的数据与真实数据进行比较，通过定量精度指标选择合适的窗口。在一定大小窗口下，利用公式(1)求得窗口内每个混合像元内各类别丰度，组成到n²行k列的丰度矩阵A(n²×k)(n²＞k)(n表示矩形窗口大小，k表示窗口内类别数)，利用公式(3)进行最小二乘法求解窗口内各类别的光谱值。然后，把求得的值依照类别对应到窗口内相应位置上，获得降尺度的数据。

公式的约束条件：

其中，式(2)和式(3)中，R(i，t)为t时期i位置低分辨率像元的反射率；f_c(i，c)为i位置低分辨率像元内类别c占该像元的面积比；为t时期类别c的平均反射率；ξ(i，t)为残差；k为窗口内类别数；

三、融合影像生成

CDSTARFM利用t₀和t_k两期低分辨率影像降尺度的数据和t_k时期中分辨率遥感影像预测t₀时期中分辨率遥感数据。STARFM算法实现的关键在于如何权重空间信息，且这个权重算法可以自适应于不同的区域。通过在一定窗口内筛选与中心像元光谱相似的相似性像元，利用相似性像元的光谱差异性，时间差异性和与中心像元的相对距离来加权计算预测期的中心像元值。由于CDSTARFM继承了STARFM算法，其实现过程也是通过加权窗口内相似性像元计算预测期的中心像元值。

窗口内中心像元的相似性像元由t_k时期中分辨率影像的green、red和NIR波段的阈值θ_b确定窗口内中心像元的相似性像元，阈值θ_b通过公式(4)计算，窗口内与中心像元差值的绝对值小于各波段θ_b的像元作为相似性像元，

式(4)中，θ_b表示窗口内b波段的阈值，N为窗口内中分辨率像元个数，x_i为i位置的像元反射率，μ为窗口内像元反射率均值，m为类别数；

相似像元的权重大小W_ijk有三个指标来衡量：1)光谱的差异性S_ijk；2)时间的差异性T_ijk；3)相似性像元与中心像元的相对距离D_ijk；

T_ijk＝|M(x_i，y_j，t_k)-M(x_i，y_j，t₀)| (6)

步骤4.3利用公式(9)加权窗口内相似像元的反射率计算预测期的中心像元反射率，生成最终的融合影像

上述公式(5)-(9)中，w为窗口大小，L(x_i，y_j，t_k)和M(x_i，y_j，t_k)分别为t_k时期给定位置(x_i，y_j)的中分辨率数据和降尺度数据的像元值，(x_w/2，y_w/2)为窗口的中心像元，光谱差异性S_ijk值越小说明给定位置与邻近像元的光谱相似度越高，赋予较高的权重；时间差异性T_ijk值越小说明该段时间内光谱变化越小，赋予较高的权重；相对距离D_ijk值越小，赋予较高的权重。

四、测试数据与预处理

测试数据覆盖的区域(37.52°-37.73°N，115.44°-115.71°E)位于河北省衡水市境内，面积约为25km×25km(50×50个MODIS像元)。实验区的Landsat8影像如图3所示，该区域以农业种植为主，作物类型以玉米为主，兼有棉花、大豆、花生等，影像右上区域为居民区，中间有一定面积的水体，地表类型复杂，在MODIS像元尺度下，纯净像元较少。

测试数据(见表1，表2)为2景2014-08-19和2014-09-04轨道号为path/row(123/034)Landsat8(OLI)影像和5景2014-06-30、07-13、07-26、08-19和09-04轨道号为h27/v05的地表日反射率产品数据(MOD09GA)。Landsat8影像的坐标系为UTM-WGS84Zone 50N，由于Landsat 8数据已进行过基于地形数据几何校正，可以直接使用。利用ENVI5.1提供的定标和FLAASH大气校正模块把DN值转换为地表真实反射率。日反射率MOD09GA数据采用HDF格式和正弦(Sinusoidal)坐标系统，由于该数据是经过大气校正和几何校正的2 级产品数据并附带有详细的地理坐标信息，可以利用MODIS重投影工具(MODIS Re-projection Tool，MRT)转换为UTM-WGS84坐标系和Geo-tif数据格式，并用最近邻域法重采样像元为480米，便于后续的混合像元分解。

表1 Landsat8和MODIS数据信息

表2 测试数据用到的Landsat8和MOD09GA对应波段的信息

五、测试数据与预处理

本发明利用2014-08-19的Landsat8多光谱数据进行ISODATA聚类，合并后类别为10类。首先，对2014-08-19的MODIS数据在窗口大小W(MODIS像元)为5×5、7×7、11×11、15×15、21×21、31×31和41×41分别进行像元分解降尺度实验，降尺度数据与真实数据通过相关系数(r)、均方根误差(RMSE)和ERGAS(the Erreur Relative GlobaleAdimensionallede Synthèse)等定量指标对各窗口下降尺度数据进行评价并选出最佳窗口，然后在最佳窗口下对2014-09-04的MODIS进行像元分解获得该期的降尺度数据。本发明提供的方法利用2014-08-19和2014-09-04两期最佳窗口下降尺度的数据与基期2014-08-19Landsat8数据进行不同窗口下green、red和NIR波段的数据融合实验，窗口w(Landsat OLI像元)设置为7×7、11×11、31×31、61×61、101×101和151×151。为了便于该算法与STARFM算法比较，利用相同的数据和窗口进行STARFM算法的数据融合实验。两种算法的融合数据与真实数据通过相关系数(r)、均方根误差(RMSE)和ERGAS进行融合数据的质量评价。

式(10)(11)中，N表示像元个数，j表示融合数据和真实数据的像元位置，x_j和y_j分别表示j位置融合数据和真实数据的像元值；和分别表示融合数据和真实数据的平均值。

式(12)中，h表示高分辨率影像的像元大小(30米)，1表示低分辨影像的像元大小(480米)；N_ban表示波段的数量；RMSE_i表示i波段融合数据与真实数据对应波段的均方根误差；M_i表示i波段真实数据的平均值。ERGAS值越小说明融合数据与真实数据在空间上的差异就越小，融合数据的质量就越高。

六、结果分析

6.1定量分析

表3、表4和表5分别表示像元分解降尺度算法、STARFM模型和本发明提供的方法在不同窗口下融合的Green、Red和NIR波段数据与同期真实Landsat8数据对应波段的相关系数(r)、均方根误差(RMSE)和ERGAS指标。

从表3中2014-08-19MODIS降尺度数据的各项评价指标可以看出，当滑动窗口W为31×31个MODIS像元时，评价指标几乎均达到了最佳值(表中的黑体数值)。因此，选择该窗口作为最佳窗口，在该窗口下对2014-09-04MODIS进行降尺度处理，其融合精度指标见表3。

表3 不同窗口2014年8月20日和9月5日MODIS数据像元降尺度数据精度指标

注：黑体数值表示实验中的最佳值

从表4中可以看出，当窗口大小为31×31个Landsat OLI像元时融合数据各波段的相关系数、均方根误差RMSE以及ERGAS指标几乎同时达到了最佳值，相关系数r达到了90％左右。本发明提供的方法融合数据的精度(表5)在11×11个Landsat8OLI像元时达到最佳。由于本发明提供的方法是利用MODIS降尺度的数据替代STARFM中重采样的MODIS数据，MODIS数据经过分解降尺度比直接重采样更能真实地反映表信息，在相同大小的窗口内，本发明提供的方法搜寻到“纯净像元”的概率增加，从而提高融合数据的质量。

与STARFM算法的最佳窗口(31×31)相比，本发明提供的方法最佳窗口(11×11)明显减小，并且也提高了融合数据的精度。说明本发明提供的方法能够保证在较小的搜索窗口内寻找到相似性像元，相同窗口大小下，本发明提供的方法的融合数据各波段的相关系数、均方根误差(RMSE)和ERGAS均高于STARFM算法。由此可见，本发明提供的方法可以在较小内融合出更高精度的数据。

图2表示三种融合方法在最佳窗口下的融合数据与真实Landsat8数据的散点图。由于像元分解降尺度算法削弱了地物的空间差异性，得到的是窗口内端元类别的平均反射率，使得散点图上会出现“条带”。而STARFM和本发明提供的方法是利用中心像元周围相似性像元的信息加权进行数据融合很好地体 现了地物的空间差异性，不会出现“条带”现象，对应波段的散点较好地分布在1∶1对角线两侧，相关系数(r)也明显高于像元分解降尺度算法。本发明提供的方法与STARFM算法相比，融合数据的相关系数(r)提高了约2％，散点图的集中性更好。由此可见，本发明提供的方法的融合数据与真实数据具有更高的相关性。

表4 STARFM算法在不同窗口w的融合数据精度指标

注：黑体数值表示实验中的最佳值

表5 本发明提供的方法不同窗口w的融合数据精度指标

注：黑体数值表示实验中的最佳值

6.2融合图像分析

像元降尺度算法、STARFM算法和本发明提供的方法在最佳窗口下融合的遥感影像分别见图3(c)-(e)。三种方法融合的影像与2014-08-19(图3(a)和2014-09-04(图3(b)真实Landsat8影像进行目视比较，图3(c)-(e)与(b)相比较可以看出，三种方法对面积较大的地物(水体、居民区和较宽的道路)均可清晰地识别。但从放大的局部细节图(h)-(j)来看，在较破碎 的区域，像元降尺度算法的融合影像(h)有“图斑”，不能识别较小的地物及地物内部的纹理信息，而(i)和(j)能够反映出细小地物和纹理特征；由于STARFM算法使用的是重采样的MODIS数据，使得重采样的像元之间具有很强的均质性，这导致融合影像中会出现MODIS像元边界(图3(i)中黑色框内所示)，而本发明提供的方法的融合影像可以有效消除MODIS边界。图3(g)中黄色椭圆框标注是(f)到(g)是发生变化的地物，三种方法融合方法均没有捕捉到该变化信息，这是因为用于融合的基期Landsat影像(图3(f))没有体现该变化信息，由于本发明提供的方法大部分继承了STARFM算法，所以也没有解决该问题；但Hilker等提出的STAARCH方法通过选择合适的基期数据解决了这一问题。总体来说，本发明提供的方法融合的影像目视效果更接近真实影像。

像元分解降尺度算法获得的是局部区域内端元的平均反射率，地物的空间差异性小，融合影像容易会出现“图斑”；STARFM算法由于在地块破碎区域很难找到纯净MODIS像元，这在一定程度上会降低融合数据的精度；而两者结合的本发明提供的方法是利用MODIS的降尺度数据替代STARFM中重采样的数据，在一定程度上更真实地反应地表真实情况，增加了纯净像元的数量，提高了数据融合的精度。

CDSTARFM比STARFM具有较高的融合精度。从表4、表5和图2可以看出，本发明提供的方法融合的green、red和NIR数据的相关系数r分别达到了0.91、0.92和0.96高于STARFM算法的0.89、0.90和0.95和降尺度算法的0.82、0.86和0.90以及更为集中的散点图；同时，本发明提供的方法最佳窗口w(11×11个OLI像元)小于STARFM算法最佳窗口w(31×31个OLI像元)，说明本发明提供的方法在较小的窗口内就寻找到了“纯净”像元，提高了融合数据的精度。

本发明提供的方法融合的影像目视效果与真实影像最接近。从图3可以看 出：三种方法融合的影像整体上差别不明显(图3)，但在破碎区域局部放大图上有很大的细节差异(图3(h)-(j))，像元分解降尺度的融合影像具有较为明显的“图斑”；STARFM算法的融合影像在地物分界处会出现MODIS像元边界的现象，而本发明提供的方法的融合数据不具有“图斑”和“MODIS边界”与真实影像最接近。

以上所述，仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可实施范围，举凡本领域熟练技术人员在未脱离本发明所指示的精神与原理下所完成的一切等效改变或修饰，仍应由本发明权利要求的范围所覆盖。

Claims (4)

1.一种构建高时空遥感数据的方法，其特征在于所述方法主要包括以下步骤：

步骤一、选取一个基期t_k的中、低分辨率遥感影像和时序低分辨率影像；

步骤二、对t_k时期中分辨率影像进行聚类，得到分类图像后进行类别丰度提取，得到类别丰度图；

步骤三、通过像元分解模型对步骤二得到的类别丰度图以及t₀和t_k两期低分辨率影像进行像元分解，得到t₀和t_k两期的降尺度影像；

步骤四、利用t_k时期中分辨率影像筛选相似像元，得到窗口内中心像元的相似性像元；

步骤五、利用步骤三得到的t₀和t_k两期低分辨率影像降尺度的影像数据和步骤四得到的相似性像元，通过加权窗口内相似性像元计算预测期的中心像元值，从而得到预测期t₀的中分辨率遥感影像；

所述步骤四和步骤五通过在一定窗口内筛选与中心像元光谱相似的相似性像元，利用相似性像元的光谱差异性，时间差异性和与中心像元的相对距离来加权计算预测期的中心像元值；

所述步骤四和步骤五具体包括：

步骤4.1由t_k时期中分辨率影像的green、red和NIR波段的阈值θ_b确定窗口内中心像元的相似性像元，阈值θ_b通过公式(4)计算，窗口内与中心像元差值的绝对值小于各波段θ_b的像元作为相似性像元；

式(4)中，θ_b表示窗口内b波段的阈值，N为窗口内中分辨率像元个数，x_i为i位置的像元反射率，μ为窗口内像元反射率均值，m为类别数；

步骤4.2计算相似像元的权重大小W_ijk，权重大小W_ijk有三个指标来衡量：光谱的差异性S_ijk，时间的差异性T_ijk和相似性像元与中心像元的相对距离D_ijk；

S_ijk＝|L(x_i，y_j，t_k)-M(x_i，y_j，t_k)| (5)

T_ijk＝|M(x_i，y_j，t_k)-M（x_i，y_j，t₀)| (6)

步骤4.3利用公式(9)加权窗口内相似像元的反射率计算预测期的中心像元反射率，生成最终的融合影像

上述公式中，w为窗口大小，L(x_i，y_j，t_k)和M(x_i，y_j，t_k)分别为t_k时期给定位置(x_i，y_j)的中分辨率数据和降尺度数据的像元值，(x_w/2，y_w/2)为窗口的中心像元，光谱差异性S_ijk值越小说明给定位置与邻近像元的光谱相似度越高，赋予较高的权重；时间差异性T_ijk值越小说明该段时间内光谱变化越小，赋予较高的权重；相对距离D_ijk值越小，赋予较高的权重。

2.根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤一中的一个基期t_k包括该时期的一景中分辨率遥感影像和一景低分辨率遥感影像。

3.根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤二中的丰度提取具体为通过公式(1)提取每个混合像元内各类别的丰度：

式(1)中，f_c(i，c)表示i位置低分辨率像元内c类端元的丰度，Q表示低分辨率像元内c端元的像元数，S表示低分辨率像元内所有端元的像元数。

4.根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤三具体包括：

步骤3.1混合像元的光谱值等于该像元内各地类的光谱值与其丰度的线性组合，可利用公式(2)表示：

步骤3.2选择合适窗口进行混合像元分解，利用公式(3)求解窗口内各类别的光谱值：

公式的约束条件：

其中，式(2)和式(3)中，R(i，t)为t时期i位置低分辨率像元的反射率；f_c(i，c)为i位置低分辨率像元内类别c占该像元的面积比；为t时期类别c的平均反射率；ξ(i，t)为残差；k为窗口内类别数；

步骤3.3把求得的各类别光谱值依照类别对应到窗口内相应像元上，获得降尺度的数据。