提供的方法融合数据的精度（表5)在11 X 11个Landsat80LI像元时达 到最佳。由于本发明提供的方法是利用MODIS降尺度的数据替代STARFM中重采样的MODIS 数据，MODIS数据经过分解降尺度比直接重采样更能真实地反映表信息，在相同大小的窗口 内，本发明提供的方法搜寻到"纯净像元"的概率增加，从而提高融合数据的质量。
[0089] 与STARFM算法的最佳窗口（31X31)相比，本发明提供的方法最佳窗口（11X11) 明显减小，并且也提高了融合数据的精度。说明本发明提供的方法能够保证在较小的搜索 窗口内寻找到相似性像元，相同窗口大小下，本发明提供的方法的融合数据各波段的相关 系数、均方根误差（RMSE)和ERGAS均高于STARFM算法。由此可见，本发明提供的方法可以 在较小内融合出更高精度的数据。
[0090] 图2表示三种融合方法在最佳窗口下的融合数据与真实LandsatS数据的散点图。 由于像元分解降尺度算法削弱了地物的空间差异性，得到的是窗口内端元类别的平均反射 率，使得散点图上会出现"条带"。而STARFM和本发明提供的方法是利用中心像元周围相似 性像元的信息加权进行数据融合很好地体现了地物的空间差异性，不会出现"条带"现象， 对应波段的散点较好地分布在I : 1对角线两侧，相关系数（r)也明显高于像元分解降尺 度算法。本发明提供的方法与STARFM算法相比，融合数据的相关系数（r)提高了约2%，散 点图的集中性更好。由此可见，本发明提供的方法的融合数据与真实数据具有更高的相关 性。
[0091 ] 表4 STARFM算法在不同窗口 w的融合数据精度指标
[0093] 注：黑体数值表示实验中的最佳值
[0094] 表5本发明提供的方法不同窗口 w的融合数据精度指标
[0096] 注：黑体数值表示实验中的最佳值
[0097] 6. 2融合图像分析
[0098] 像元降尺度算法、STARFM算法和本发明提供的方法在最佳窗口下融合的遥感影像 分别见图3 (c) - (e)。三种方法融合的影像与2014-08-19 (图3 (a)和2014-09-04 (图3 (b) 真实LandsatS影像进行目视比较，图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方法通过选择合适的基期数据解决了这一问 题。总体来说，本发明提供的方法融合的影像目视效果更接近真实影像。
[0099] 像元分解降尺度算法获得的是局部区域内端元的平均反射率，地物的空间差异 性小，融合影像容易会出现"图斑"；STARFM算法由于在地块破碎区域很难找到纯净MODIS 像元，这在一定程度上会降低融合数据的精度；而两者结合的本发明提供的方法是利用 MODIS的降尺度数据替代STARFM中重采样的数据，在一定程度上更真实地反应地表真实情 况，增加了纯净像元的数量，提高了数据融合的精度。
[0100] ⑶STARFM比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(llX 11个OLI像元）小于STARFM算法最佳窗 口 w(31X31个OLI像元），说明本发明提供的方法在较小的窗口内就寻找到了"纯净"像 元，提高了融合数据的精度。
[0101] 本发明提供的方法融合的影像目视效果与真实影像最接近。从图3可以看出：三 种方法融合的影像整体上差别不明显（图3)，但在破碎区域局部放大图上有很大的细节差 异（图3 (h) - (j))，像元分解降尺度的融合影像具有较为明显的"图斑";STARFM算法的融合 影像在地物分界处会出现MODIS像元边界的现象，而本发明提供的方法的融合数据不具有 "图斑"和"MODIS边界"与真实影像最接近。
[0102] 以上所述，仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可 实施范围，举凡本领域熟练技术人员在未脱离本发明所指示的精神与原理下所完成的一切 等效改变或修饰，仍应由本发明权利要求的范围所覆盖。
【主权项】
1. 一种构建高时空遥感数据的方法，其特征在于所述方法主要包括以下步骤： 步骤一、选取一个基期（tk)的中、低分辨率遥感影像和时序低分辨率影像； 步骤二、对tk时期中分辨率影像进行聚类，得到分类图像后进行类别丰度提取，得到类 别丰度图； 步骤三、通过像元分解模型对步骤二得到的类别丰度图以及〖。和tk两期低分辨率影像 进行像元分解，得到1。和tk两期的降尺度影像； 步骤四、利用tk时期中分辨率影像筛选相似像元，得到窗口内中心像元的相似性像 元； 步骤五、利用步骤三得到的t。和、两期低分辨率影像降尺度的影像数据和步骤四得到 的相似性像元，通过加权窗口内相似性像元计算预测期的中心像元值，从而得到预测期t。 的中分辨率遥感影像。2. 根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤一中的 一个基期（tk)包括该时期的一景中分辨率遥感影像和一景低分辨率遥感影像。3. 根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤二中的 丰度提取具体为通过公式（1)提取每个混合像元内各类别的丰度：式⑴中，fji，c)表示i位置低分辨率像元内c类端元的丰度，Q表示低分辨率像元 内c端元的像元数，S表示低分辨率像元内所有端元的像元数。4. 根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤三具体 包括： 步骤3. 1混合像元的光谱值等于该像元内各地类的光谱值与其丰度的线性组合，可利 用公式⑵表示。其中，式⑵和式（3)中，R(i，t)为t时期i位置低分辨率像元的反射率；fji，c)为i位置低分辨率像元内类别c占该像元的面积比；为t时期类别c的平均反射率； Ui，t)为残差；k为窗口内类别数； 步骤3. 3把求得的各类别光谱值依照类别对应到窗口内相应像元上，获得降尺度的数 据。5. 根据权利要求1所述的构建高时空遥感数据的方法，其特征在于：所述步骤四和步 骤五通过在一定窗口内筛选与中心像元光谱相似的相似性像元，利用相似性像元的光谱差 异性，时间差异性和与中心像元的相对距离来加权计算预测期的中心像元值。6. 根据权利要求5所述的构建高时空遥感数据的方法，其特征在于：所述步骤四和步 骤五具体包括： 步骤4. 1由巩时期中分辨率影像的green、red和NIR波段的阈值0 !^角定窗口内中心 像元的相似性像元，阈值叭通过公式（4)计算，窗口内与中心像元差值的绝对值小于各波 段9b的像元作为相似性像元。式（4)中，9b表示窗口内b波段的阈值，N为窗口内中分辨率像元个数，Xi位置的 像元反射率，y为窗口内像元反射率均值，m为类别数； 步骤4. 2计算相似像元的权重大小Wljk，权重大小Wljk有三个指标来衡量：光谱的差异 性Sljk，时间的差异性Tljk和相似性像元与中心像元的相对距离Dljk;步骤4. 3利用公式（9)加权窗口内相似像元的反射率计算预测期的中心像元反射率， 生成最终的融合影像上述公式中，w为窗口大小，LU;,y_j，tk)和M(Xi，y_j，tk)分别为tk时期给定位置（X;，y])的中分辨率数据和降尺度数据的像元值，（xw/2，yw/2)为窗口的中心像元，光谱差异性S1]k 值越小说明给定位置与邻近像元的光谱相似度越高，赋予较高的权重；时间差异性T1 ]k值 越小说明该段时间内光谱变化越小，赋予较高的权重；相对距离D1 ]k值越小，赋予较高的权 重。
【专利摘要】本发明公开了一种构建高时空遥感数据的方法，结合像元分解降尺度方法和STARFM模型的CDSTARFM方法，首先利用像元分解降尺度算法对低分辨率数据进行降尺度处理，然后用降尺度数据取代STARFM算法中直接重采样的低分辨率数据，最后利用两者相结合的CDSTARFM算法进行数据融合。该方法可以有效解决像元分解降尺度融合数据的“图斑”现象和STARFM模型寻找纯净MODIS像元难的问题。利用本发明提供的方法融合Landsat8和MODIS数据构建高时空分辨率遥感数据。结果表明，本发明提供的方法比STARFM和像元分解降尺度算法具有更高的融合精度。
【IPC分类】G06T3/40, G06T5/50
【公开号】CN105046648
【申请号】CN201510354552
【发明人】张锦水, 谢登峰, 潘耀忠, 袁周米琪, 云雅, 孙佩军
【申请人】北京师范大学
【公开日】2015年11月11日
【申请日】2015年6月25日