CN109840539A - 一种基于地块图斑的遥感时空数据融合方法 - Google Patents

Abstract

本发明公布了一种基于地块图斑的遥感时空数据融合方法，在获取地表真实地块图斑基础上，以图斑为约束条件，利用研究区的一期高分影像以及多期中分辨率影像，通过对高分影像进行无监督分类获取地表覆盖类型并且计算每个图斑的反射率；以图斑为约束条件对中分辨率数据进行最优解混；在对中分辨率影像进行解混的基础上，推测中分辨率数据对应的每期高分辨数据。本发明可以提高传统的基于像元级别的时空数据融合方法，与采用学习模型融合方法相比，可以降低难度，提高稳定性。

Description

一种基于地块图斑的遥感时空数据融合方法

技术领域

本发明提出了一种以地块图斑为基本处理单元，并以地块图斑边界作为约束知识，通过高分影像与中分时间序列影像的协同处理实现地块尺度时空数据融合的方法。

背景技术

遥感影像的空间分辨率和时间分辨率是一对矛盾体：低空间分辨率影像由于卫星轨道较高，因而覆盖面宽，重访周期短（卫星经过同一地面点），时间分辨率高；高空间分辨率影像由于卫星轨道低，因而覆盖范围较窄，重访周期长，时间分辨率低。单一卫星数据源难以同时满足高空间和高时间分辨率的应用需求。

随着遥感技术的发展，遥感观测数据日益丰富，数据类型覆盖了高、中、低多种空间分辨率，遥感应用的范围逐渐扩大，同时也朝着精细化方向不断发展。精细化的遥感应用反过来对高空间、高时间分辨率数据提出了更高要求。遥感时空数据融合技术为解决遥感数据时空分辨率的矛盾性提供了有效的解决途径。时空数据融合通过将高空间-低时间分辨率与低空间-高时间分辨率两种卫星影像数据的协同处理，生成同时具有高空间分辨率和高时间分辨率的合成影像。

在方法上，时空数据融合目前主要有3种方式：（1）基于数学变换模型的融合。该方法首先对二维影像数据进行相应的数学变换，然后对变换后的数据进行处理，最后进行反变换获得期望影像，一般采用的变换方法为小波变换。（2）基于重建模型的融合。该方法采用同一区域的1-2期高空间-低时间分辨率影像与该时段内的多幅低空间-高时间分辨率的时间序列影像（其中首尾2期与高分数据的时相相同）进行融合。首先根据同一时相的高分和低分影像计算不同分辨率像元间的对应关系，获得高-低分辨率之间的像元换算函数f(r)；然后对低空间-高时间分辨率的时间序列影像像元进行分解并计算像元反射率随时间的变化函数关系f(t);最后，将低分像元时间变化关系f(t)与高-低分辨率像元换算函数f(r)耦合，进而插值得到时间序列低分影像对应时相的高空间分辨率合成影像。STARFM时空自适应融合是这类方法的典型代表，也是应用最多的一类算法，其实质是对影像的光谱信息和空间信息进行插值，融合结果与真实结果仍存在一定的差距。（3）基于学习模型的融合。该方式通常采用压缩感知和稀疏表达的方法对遥感影像进行预处理，在此基础上基于机器学习的手段进行训练与拟合。该方法通常需要大量时间用于模型训练，受所选用的样本影响，模型鲁棒性和可拓展性有待提高，并且对不同类型影像的分辨率差值有要求，尚不适用于大区域的遥感数据时空融合。相关文献：刘建波, 马勇, 武易天, 等. 遥感高时空融合方 法的研究进展及应用现状[J]. 遥感学报, 2016, 20(5): 1038-1049.Malenovský Z, Bartholomeus H M, Acerbi-Junior F W, et al. Scaling dimensions in spectroscopy of soil and vegetation[J]. International Journal of Applied Earth Observation and Geoinformation, 2007, 9(2): 137-164. Gao F, Masek J, Schwaller M, et al. On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance[J]. IEEE Transactions on Geoscience and Remote sensing, 2006, 44(8): 2207-2218. Huang B, Song H. Spatiotemporal reflectance fusion via sparse representation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(10): 3707- 3716.Zeyde R, Elad M, Protter M. On single image scale-up using sparse- representations[C]//International conference on curves and surfaces. Springer, Berlin, Heidelberg, 2010: 711-730.

上述时空数据融合方法均为在像素尺度上进行，主要针对低空间分辨率影像（主要为MODIS数据，250m或500m分辨率，重访周期8天）和中空间分辨率影像（主要为LANDSAT数据，30m分辨率，重访周期16天），通过不同的融合模型获得每8天的30m分辨率合成影像。而针对高空间分辨率（米级／亚米级）和多种来源的中空间分辨率数据的融合方法未见报道。随着遥感应用朝着精细化方向发展，尤其是精准农业的应用需求，地块级的高空间-高时间分辨率观测信息成为了限制应用发展的瓶颈。

本发明提供了一种基于地块图斑的遥感时空数据融合方法，在米级／亚米级地块边界数据支持下，实现高空间分辨率和多种来源的中空间分辨率数据的融合，融合结果以地块对象的反射率特征来表达，同时这种以地块对象为基本处理单元的融合方法更为贴近应用需求，可以更为便捷地融入外部知识和环境信息，用于提高后继应用分析的几何与属性精度。相比于前述传统的融合方式，本方法提高了融合精度，同时具有更好的稳定性和可操作性。

发明内容

本发明方法的输入数据包括地块图斑矢量、1-2期高空间分辨率影像（米级或亚米级）和中等空间分辨率的时间序列影像（16-250米均可，可来自不同卫星数据）。方法的前提是输入的影像数据已经过几何校正及辐射校正，并假设融合时段内影像对应区域的地物类型不发生变化（如耕地变为建筑物）以及相同地类的时间变化率相同。

本发明方法的理论依据是线性混合像元解混方法，主要思路是以地块边界为约束条件，基于高分影像获取高精度的地类分布信息与t ₀时刻地块图斑的真实反射率r ₀；在此基础上，结合地类分布信息，对时间序列(t ₀、t ₁…t _n)中分影像进行线性光谱解混获得各地类的时序平均反射率r ₀’、r ₁’…r _n’，进一步计算获得不同时相间的反射率增量 QUOTE 、 QUOTE … QUOTE ；进一步，结合高分影像获得的地块图斑类型与t ₀时刻地块图斑的真实反射率r ₀，计算获得各个时相的高分辨率地块图斑反射率预测值（融合结果）r ₁、r ₂…r _n；进一步可计算图斑的多种光谱特征时序，同时基于高分影像还可获得图斑对象的多种纹理特征。

数学推导过程如下：已知同区域时刻的高分和中分反射率影像对，地块图斑矢量数据以及时刻的中分反射率影像，求取时刻各个地块图斑内包含的高分像元平均反射率（即 地块图斑反射率）。

根据混合像元线性分解理论，中分辨率像元的反射率等于该像元内各类地物的反射率(高分辨率像元反射率)与其面积占比的加权和：

（1）

其中： 为第i个中分像元在b波段的发射率，为c类地物在b波段的平均反射率，为c类地物占中分像元的百分比（丰度），e(i)为误差项。若已知和各类地物的面积占比，则各类地物的平均反射率可采用最小二乘法进行方程解算。为了减少地类光谱的空间变异性影响，一般情况下在中分影像上采用一定大小的窗口进行方程求解，此时各地类的可从高分辨率的影像分类图中计算获取。

在给定大小窗口内，地类c从到时刻的平均反射率的变化量可表示为公式（2），其中p为窗口内c类地物高分像元的个数。

（2）

在同一地类的时间变化率相同的前提假设下，由公式（2）可得到式（3），即地类c从到时刻的平均反射率变化量等于同时段内属于c类地物高分像元的反射率变化量。

（3）

由式（3）可得下式（4），表明时刻高分像元反射率可由时刻高分像元反射率与地类c到时刻的平均反射率变化量求和得到，而地类c的平均反射率变化量可根据公式（1）对和时刻的中分影像进行基于最小二乘法的混合像元分解求取。

（4）

以上是像素级时空数据融合的推导过程。对于地块图斑而言，其地块内高分像元为同一地物类型，则地块图斑反射率为：

（5）

式中q为地块图斑所包含的高分像元数目。从到时刻的地块图斑反射率变化量可表示为：

（6）

在同一地类的时间变化率相同的前提假设下，由公式（6）可得到式（7）

（7）

式（7）等价于下式：

（8）

式（8）表明，时刻的地块图斑反射率等于时刻的地块图斑反射率与地类c从到时刻的平均反射率变化量之和。其中时刻的地块图斑反射率可从高分影像得到，而地类c的平均反射率变化量可从和时刻的中分影像获得。

以上是基于地块图斑的时空数据融合推导过程。在实际计算过程中，

具体步骤如下：

1） 以地块图斑作为约束条件对高分影像进行最优类别估计的非监督分类，获取地类分布信息和地块图斑所属的类型，同时取地块图斑内部像元的反射率均值作为t ₀时刻的图斑真实反射率。

2） 计算中分影像像元内部的地类面积比率，以图斑为约束条件，采用迭代逼近的方式获得局部光谱解混的最佳窗口，实现窗口内各地类平均反射率的最优分解对，实现中分影像的进行解混。

3） 利用第2）以及第1）的结果，推测中分数据对应的每期高分数据的反射率增量，在此基础上获取各期具有高分特征的数据。

上述步骤特征在于：

步骤1）以图斑为约束条件对高分数据进行无监督分类，获取地表覆盖类型及图斑级别的反射率。

步骤2）以图斑为约束条件对中分数据进行解混，获得每期高分数据反射率的增量。

步骤3）实现图斑级别的高分数据推理，获取对应中分影像的高分时间序列影像。

本发明与以往的中分—高分数据时空数据融合方法相比具有以下特点：以图斑为约束条件对高分数据进行非监督分类获取地表覆盖类型的过程；以图斑为约束条件求取图斑级别的地表真实反射率；以图斑为约束条件对中分影像进行解混的过程；以图斑为约束条件实现对应每期中分影像的高分影像推理；以图斑为约束条件的时间序列植被指数计算及影像纹理提取。

附图说明

图1是基于地块图斑的遥感时空数据融合方法示意图；

图2是低分时序影像中得到高精度的时间变化信息示意图；

图3是分类状态示意图，A、B均为过分类，C为合适分类；

图4是高分影像最优分类示意图；

图5是中分影像像元解混示意图，设地块半径为R，像元半径为 r,像元与地块中心的距离为d，当满足r+d<R时，则本地块不参与中分影像混合像元分解计算，提高了像元解混的速度；

图6是基于图斑的时空融合方法与传统方法的对比示意图，图中圈出来位置为本方法融合后数据与传统方法融合后数据分别进行分类的结果对比；

图7是基于图斑的时空融合方法与传统方法提取植被指数NDVI精度对比示意图，图中横坐标为影像获取的NDVI，纵坐标为拟合NDVI。

具体实施方式

图1 示意了本发明的主要实现思路。首先对获取的研究区的中分和高分影像进行几何校正以及辐射校正，在此基础上以图斑为约束条件，进行高分数据的地表覆盖类划分以及图斑级别的地表真实反射率计算，并以此进行中分数据的解混，进一步的实现中分数据对应的每期的预测高分数据；最后通过预测的高分数据进行时间序列植被指数计算。应用本发明方法对广西崇左市江州区进行了中-高分时空数据的融合试验。

图2示意了由低分时序影像进高分影像预测过程。A1，A2，A3为江洲区三期中分影像，B1为获取的一期高分影像，缺失B2，B3两期高分影像。B2’，B3’为通过A1，A2，A3三期中分影像和B1一期高分影像进行预测获得的对应的B2，B3高分影像。其算法是：

H(x _i ,y _j ,t _k )= QUOTE

式中H(x _i ,y _j ,t _k )为预测高分像元值，H(x _i ,y _j ,t ₀ )为t ₀ 时期真实高分像元值，L(x _i ,y _j , t _k )和L(x _i ,y _i ,t ₀ )分别为t _k 和t ₀ 时期的中分像元值x _i ，，w为搜索窗口，W _ijk 为联合权重矩阵。

图4为高分影像最优分类示意图。该分类结果借助于已有的参考图斑，通过对高分影像进行迭代得到，其最优判别依据如下：

SCR= QUOTE

SCAR= QUOTE

当SCR>0.9&&SCAR>0.95时，确认本发明达到最优分类。

经过融合，试验结果如图6所示，图6A为面向对象的融合方法后地块提取结果，B为基于本发明的算法融合后进行的地块提取。经检验，相对传统方法，本发明精度随分辨率和地块大小变化稳定，相对误差在4%以下。

图 7为本发明方法融合结果与传统方法融合结果进行植被指数NDVI后的精度对比，传统方法随分辨率增大，地块读值与真实值偏差也越大；而本文方法随分辨率增加表现较为稳定，且集中于1:1线上，相对于经典方法，本方法的反射率在绝对精度上可以最高提高6—8倍。

Claims (3)

1.一种基于地块图斑的遥感时空数据融合方法，其特征在于包括以下几个步骤：

1）在高分辨率遥感数据基础上,以地块图斑为约束条件进行无监督分类，获取地表覆盖类型，并计算图斑内部的反射率；

2）以地块图斑为约束条件，对中空间分辨率数据（简称中分数据）采用迭代逼近获取局部光谱解混的最佳窗口，该最佳窗口的确定由地块所属地类的局部空间分布特征决定，然后进行像元解混，解混是指将像元分解为不同端元，并求得每一端元的丰度；

3）利用步骤1）与步骤2）的结果推测每期中分影像对应的高分数据的反射率增量，以获取的高分数据为基础，利用每个增量计算对应的具有高分特征的数据。

2.如权利要求1所述一种基于地块图斑的遥感时空数据融合方法，其特征是以图斑为约束条件对高分影像进行无监督分类获取地表覆盖类型和图斑内部平均反射率作为本图斑的反射率。

3.如权利要求1所述一种基于地块图斑的遥感时空数据融合方法，其特征是以图斑为约束条件通过计算，获取每期中分影像对应的高分数据的反射率增量，从而获得每期具有高空间分辨率特征的数据。