一种构建高时空遥感数据的方法
【技术领域】
[0001] 本发明涉及一种构建高时空遥感数据的方法，尤其涉及一种在地块破碎区域构建 高时空遥感数据的方法，属于遥感影像数据处理领域。
【背景技术】
[0002] 高空间、时间分辨率的遥感数据在监测土地覆盖变化和作物生长以及识别地物类 型等方面具有重要的作用。然而，目前获取到的卫星传感器数据存在空间分辨率、时间分辨 率、光谱分辨率相互矛盾，即同时满足高空间、高时间和高光谱分辨率的卫星遥感数据是不 现实的。比如，Landsat卫星的多光谱影像空间分辨率是30米，其适中的空间分辨率和易 获取性在植被指数提取、监测土地覆盖动态变化和生态系统方面具有广泛的应用，但因其 16天的重访周期和阴雨天气的影响，难以保证获得连续有效的影像进行农作物生长、物候 变化和自然灾害的监测；搭载在Terra/Aqua卫星上的中分辨率成像光谱仪（MODIS)能提供 时间分辨率为1天、空间分辨率为250-1000米的遥感影像，被广泛应用于大尺度区域的地 表覆盖类型监测，但MODIS数据的低空间分辨率限制其在景观破碎和异质性较强区域的应 用。若能综合中高分辨率影像（如Landsat 0LI)的高分辨率特性和MODIS时间分辨率的 优势，将会大幅度提高植被遥感监测能力。
[0003] 近些年，国内外学者提出了系列获取高时空分辨率数据的遥感数据融合模型。这 些融合模型总体分为两大类，一类是像元分解降尺度算法（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相比增加 了数据获取的难度。
【发明内容】
[0004] 本发明针对STARFM模型因在破碎区域难搜寻到纯净像元而降低融合数据精度的 问题，首先利用像元分解降尺度算法对低分辨率数据进行降尺度处理，然后将降尺度数据 取代STARFM算法中直接重采样的低分辨率数据，最后再利用两者相结合的CDSTARFM算法 进行数据融合。
[0005] 具体地，本发明提出了一种结合像元分解降尺度算法和STARFM模型的 CDSTARFM(Combining Downscale Mixed Pixel with Spatial and Temporal Adaptive Reflectance Fusion Model)算法并用于融合Landsat8和MODIS反射率数据构建每天的 Landsat8反射率数据实验。
[0006] 本发明提供一种构建高时空遥感数据的方法，主要包括以下步骤：
[0007] 步骤一、选取一个基期（tk)的中、低分辨率遥感影像和时序低分辨率影像；
[0008] 步骤二、对tk时期中分辨率影像进行聚类，得到分类图像后进行类别丰度提取，得 到类别丰度图；
[0009] 步骤三、通过像元分解模型利用步骤二得到的类别丰度图以及丨。和、两期低分辨 率影像进行像元分解，得到1。和t k两期的降尺度影像；
[0010] 步骤四、利用tk时期中分辨率影像筛选相似像元，得到窗口内中心像元的相似性 像元；
[0011] 步骤五、利用步骤三得到的〖。和、两期低分辨率影像降尺度的影像数据和步骤四 得到的相似性像元，通过加权窗口内相似性像元计算预测时期的中心像元值，从而得到预 测期t。的中分辨率遥感影像。
[0012] 优选的，上述步骤一中的一个基期（tk)包括该时期的一景中分辨率遥感影像和一 景低分辨率遥感影像。
[0013] 优选的，上述步骤二中的丰度提取具体为通过公式（1)提取每个混合像元内各类 别的丰度：
[0015] 式（1)中，fji，C)表示i位置低分辨率像元内C类端元的丰度，Q表示低分辨率 像元内c端元的像元数，S表示低分辨率像元内所有端元的像元数。
[0016] 优选的，上述步骤三具体包括：
[0017] 步骤3. 1混合像元的光谱值等于该像元内各地类的光谱值与其丰度的线性组合， 可利用公式（2)表示。
[0019] 步骤3. 2选择合适窗口进行混合像元分解，利用公式（3)求解窗口内各类别的光 谱值。
[0022] 其中，式⑵和式⑶中，R(i，t)为t时期i位置低分辨率像元的反射率；fji， c)为i位置低分辨率像元内类别c占该像元的面积比为t时期类别c的平均反射 率；ξ (i，t)为残差；k为窗口内类别数；
[0023] 步骤3. 3把求得的各类别光谱值依照类别对应到窗口内相应地类的像元上，获得 降尺度的数据。
[0024] 优选的，上述步骤四和步骤五通过在一定窗口内筛选与中心像元光谱相似的相似 性像元，利用相似性像元的光谱差异性，时间差异性和与中心像元的相对距离来加权计算 预测期的中心像元值。
[0025] 优选的，上述步骤四和步骤五具体包括：
[0026] 步骤4. 1由tk时期中分辨率影像的green、red和NIR波段的阈值Θ b确定窗口内 中心像元的相似性像元，阈值θb通过公式（4)计算，窗口内与中心像元差值的绝对值小于 各波段9b的像元作为相似性像元。
[0028] 式（4)中，0b表示窗口内b波段的阈值，N为窗口内中分辨率像元个数，X 1Si 位置的像元反射率，μ为窗口内像元反射率均值，m为类别数；
[0029] 步骤4. 2计算相似像元的权重大小Wljk，权重大小Wljk有三个指标来衡量：光谱的 差异性Sljk，时间的差异性Tljk和相似性像元与中心像元的相对距离D ljk;
[0034] 步骤4. 3利用公式（9)加权窗口内相似像元的反射率计算预测期的中心像元反射 率，生成最终的融合影像
[0035]
[0036] 上述公式中，w为窗口大小，L(Xi，y_j，tk)和Μ(Χρ y_j，tk)分别为tk时期给定位置 (Xl，y])的中分辨率数据和降尺度数据的像元值，（xw/2，y w/2)为窗口的中心像元，光谱差异性 S1]k值越小说明给定位置与邻近像元的光谱相似度越高，赋予较高的权重；时间差异性T 1]k值越小说明该段时间内光谱变化越小，赋予较高的权重；相对距离D1]k值越小，赋予较高的 权重。
[0037] 本发明提供的构建高时空遥感数据的方法可以有效解决像元分解降尺度融合 数据的"图斑"现象和STARFM模型寻找纯净MODIS像元难的问题。利用本发明提供的 方法融合Landsat8和MODIS数据构建高时空分辨率遥感数据。结果表明，本发明提供 的方法比STARFM和像元分解降尺度算法具有更高的融合精度及更接近真实影像的目视 效果，如NIR波段，相关系数r(0·96vs0·95vs0·90;RMSE(0·0245vs0·0300vs0·0401) ;ERGAS(0. 5416vs 0.6507vs 0.8737)〇
【附图说明】
[0038] 图1为本发明流程示意图；
[0039] 图2为三种不同方法预测的反射率与真实反射率的散点图；
[0040] 图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)分别表示对应的子区域影像。
【具体实施方式】
[0041] 为了便于本领域普通技术人员理解和实施本发明，下面结合附图及【具体实施方式】 对本发明作进一步的详细描述。
[0042] 本发明结合像元分解降尺度算法和STARFM模型的⑶STARFM(Combining Downscale Mixed Pixel with Spatial and Temp