一种构建高空间分辨率时序遥感数据的方法

1.本发明涉及遥感影像处理技术领域，尤其涉及一种构建高空间分辨率时序遥感数据的方法。


背景技术：

2.遥感技术以其空间覆盖范围广、时间周期固定等优势为区域地表信息动态监测提供了重要的数据来源。然而，受到传感器本身和各种不利天气条件的限制，目前单一传感器很难兼具空间和时间分辨率。高时间分辨率的影像其空间分辨率较低，混合像元的影响制约了地表信息监测的精度；而高空间分辨率的影像回访周期长、时效性差，难以满足连续的动态监测。遥感数据时空融合技术是满足高时空应用需求的有效方式。遥感影像融合是指利用一对或多对准同步的粗-细分辨率影像作为基准期影像，结合待预测日期的粗分辨率影像估测待预测日期的细分辨率影像。
3.针对特定的应用场景学者们提出了诸多的时空融合算法，其中fit-f方法仅仅需要起始日期的一对粗-细分辨率的基准期影像，并且计算效率高，对于有物候变化的场景表现出较好的预测精度，因此适用于地表信息连续的动态监测。但是该方法存在两方面的不足：一方面，受粗分辨率影像上混合像元的影响，该方法在异质性强但变化强度低的区域表现不佳，不擅长捕捉影像的结构和纹理；另一方面，由于仅利用一对基准期影像，对于基准期和预说明书测期的细分辨率影像间存在变化，而粗分辨率影像上不可见的情况，不能被捕捉到。因此，有必要改进现有的fit-fc方法，使其更加适应于连续的地表信息动态监测。
4.此外，sentinel-2a/b卫星具有5天的重访周期，即使受云等天气条件的影响造成数据缺失，其仍然为地表信息监测提供了较为密集的高空间分辨率影像。如何充分利用可获得的多幅高空间分辨率影像提供的信息，构建出更高质量的高空间分辨率时序遥感数据，也是亟待解决的问题。


技术实现要素：

5.为了克服上述现有技术的局限，本发明提供了一种构建高空间分辨率时序遥感数据的方法。
6.本发明的技术方案如下：一种构建高空间分辨率时序遥感数据的方法，具体包含以下步骤：s1、低空间分辨率影像的混合像元分解：基于高空间分辨率影像的分类图，获取低空间分辨率影像上逐像元的各类地物的丰度，并利用像元分解模型计算每类地物的光谱反射率，获得粗影像的降尺度影像；s2、基准期与预测期的回归模型拟合：基于步骤s1中得到的降尺度影像，建立基准
期和预测期（）间反射率的线性回归模型，并应用于细影像，得到细影像的反射率拟合结果；说明书s3、残差的分配和补偿：将上述步骤s2得到的细影像的反射率升尺度至低空间分辨率，计算其与粗影像反射率的差值，即残差，并采用空间插值方法将残差分配至高空间分辨率，并补偿于的细影像；s4、拟合结果的空间滤波：在一定窗口内筛选中心像元的光谱相似像元，利用相似性像元的光谱差异性和与中心像元的相对距离来加权计算预测期的中心像元值，从而得到空间滤波后的的细影像；s5、正反向估测结果的加权：设置起始日期（）和结束日期（）两个基准期，将上述步骤s1-s4分正、反两个方向进行，采用线性加权方法综合正、反向估测结果作为加权时序融合结果；s6、结合高空间分辨率影像的连续校正：基于上述步骤s5获得的时序影像融合结果，构建地表信息的动态变化模型；利用可获得的多幅已知细影像对动态模型进行修正，使动态模型模拟值与已知细影像值的误差最小，最终得到校正后的时空融合结果。
7.优选的，步骤s1中具体包括：s11：对获得的遥感影像进行大气校正、几何校正、投影转换、重采样、质量控制等预处理；s12：采用监督分类或非监督分类法对细影像进行分类，利用公式（1）计算粗影像上每个像元内各类地物的丰度；说明书式（1）中，表示i位置处混合像元内类别q的丰度，q为混合像元中类别q的细像元数，s为混合像元中所有类别的细像元总数；s13：混合像元的反射率值可以表示为像元中各类别的反射率与其丰度的线性组合，混合像元分解模型具体为：
约束条件:式（2）中，为第i个混合像元的反射率，表示第i个混合像元中类别q的平均反射率，m为类别数，为残余误差。
8.s14：将上述步骤计算得到的各类别的平均反射率按照类别对应到相应地类的像元上，从而得到粗影像的降尺度数据；优选的，步骤s2中具体包括：根据公式（3）利用基准期和预测期的降尺度数据建立线性回归模型，需要在一定大小的移动窗口m内进行，利用最小二乘法求得移动窗口内中心像元的回归系数a和b，并应用于细影像，根据公式（4）估测得到的细影像；的细影像；说明书其中，和分别为移动窗口m内和时期降尺度数据的反射率，为拟合残余误差，和分别表示和时期细影像的反射率。
9.优选的，步骤s3中具体包括：采用平均聚合的方法将的细影像升尺度至低空间分辨率，利用公式（5）计算残差rd；将残差分配至高空间分辨率采用的空间插值方法可以选用适当的克里金插值法，并将高空间分辨率的残差rd补偿于与上述步骤s2得到的细影像的反射率；式（5）中，表示低空间分辨率下的残差，和分别为和时期粗像元的反射率，和分别为和时期粗像元中细像元的观测和估测反射率。
10.优选的，步骤s4中具体包括：根据公式（6），综合基准期细影像的各波段相似像元和中心像元的光谱差异d，在一定窗口内，差异最小的前s个像元被筛选为相似像元；根据公式（7）计算邻近像元和中心像元的空间距离，根据公式（8）确定各相似像元的权重；根据公式（9）通过加权窗口内相似性像元逐波段计算预测时期的中心像元值，从而说明书得到时期空间滤波后的细影像；时期空间滤波后的细影像；时期空间滤波后的细影像；时期空间滤波后的细影像；其中，和分别表示细影像第b波段的i位置和中心位置像元的反射率，c为用于融合的影像波段数；w为窗口半径，为高空间分辨率下的残差分配值，为w
×
w窗口内中心像元的正向预测结果。
11.优选的，步骤s5中具体包括：依据细影像的分类图，确定正反向估测结果中每一地类的权重；具体地，以正向估测过程为例，根据公式（10），由估测的细影像，确定估测影像与实际观测影像间的回归系数k和l，并应用于的正向估测影像，根据公式（11）得到的似“真实影像”；根据公式（12）和（13）计算的正向估测影像和似“真实”影像的差距，得到正向估测的相对精度；
说明书说明书其中，和分别表示的正向估测和真实观测值，和分别表示的正向估测和似“真实”观测值，表示正向估测与真实观测值之间的差距，表示正向估测的相对精度；相似地，依照上述流程，可以获得反向估测的相对精度；根据公式（14）和（15）确定各地类正反向估测的权重；根据公式（16）获得正反向估测加权的时序影像融合结果；结果；结果；。
12.优选的，步骤s6中具体包括：采用s-g滤波方法对时序融合数据进行滤波用于构建地表信息动态模型；采用一种连续校正的cacao方法（公式（17）），利用多幅已知的细影像修正动态模型，将修正后的各波段反射率作为最终的时序数据融合结果；式（17）中，bp
ture
表示已知的细影像的各波段反射率，n为可获得的细影像的数量，bp
pm
表示动态模型模拟的各波段反射率，说明书scale和shift为需要计算的修正参数。
13.本发明采用以上技术方案，与现有技术相比，具有以下技术效果：1、本发明结合了混合像元分解模型，极大地减弱了粗影像上混合像元对融合结果造成的空间模糊现象，得到的融合结果具有更丰富的纹理特征和更好的目视效果。
14.2、本发明同时考虑了基准期和预测期间拟合回归模型的残差和不同传感器间的误差，并且采用空间插值方法代替传统的插值法来分配残差，取得了更高的融合精度。
15.3、本发明综合了起始日期和结束日期对预测期的双向估测结果，有利于捕捉起始期和预测期间发生的地表覆盖类型变化，使得融合结果与实际影像更接近。
16.4、本发明充分利用了可获得的多幅已知高分辨率影像，降低了融合结果与实际影像间的误差，取得更好的时序融合结果。
附图说明
17.图1为本发明实施例的实现流程图；图2为可获得的sentinel-2影像的时间分布图；图3不同方法融合影像反射率与实际观测影像反射率的散点密度图（以2019年8月4日为例）；图4为不同方法融合影像的空间特征图（以2019年8月4日为例）。
18.说明书。
具体实施方式
19.为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明的具体实施方式作进一步详细描述。显然，以下实施例用于说明本发明，但不用来限制本发明的范围。
20.一、实施案例如图1所示，一种构建高空间分辨率时序遥感数据的方法，具体包括以下步骤：s1、低空间分辨率影像的混合像元分解；采用监督分类或非监督分类法对细影像进行分类，基于细影像的分类图，获取粗影像上逐像元的各类地物的丰度；利用像元分解模型计算每类地物的光谱反射率，获得粗影像的降尺度影像；s11、根据获取到的遥感数据级别情况，进行大气校正、几何校正、投影转换、重采样、质量控制等预处理；s12、采用监督分类或非监督分类法对基准期的细影像进行分类，根据研究区地表异质性确定分类的类别数；利用公式（1）计算粗影像上每个像元内各类地物的丰度。
21.式（1）中，表示i位置处混合像元内类别q的丰度，q为混合像元中类别q的细像元数，s为混合像元中所有类别的细像元总数；说明书s13、忽略不同传感器之间的地理配准和大气校正误差，混合像元的反射率值可以表示为像元中各类别的反射率与其丰度的线性组合，并且假设类别和丰度在基准期和预测期间不发生变化，在一定大小的窗口内，采用带约束条件的最小二乘法求解各类别的平均反射率。混合像元分解模型具体为：约束条件:
式（1）中，为第i个混合像元的反射率，表示第i个混合像元中类别q的平均反射率，m为类别数，为残余误差；s14、将上述步骤计算得到的各类别的平均反射率按照类别对应到相应地类的像元上，从而得到粗影像的降尺度数据，依次对所有预测期的粗影像进行混合像元分解。
22.s2、基准期与预测期的回归模型拟合。基于步骤1中得到的降尺度影像，建立基准期和预测期间反射率的线性回归模型，并应用于细影像，得到细影像的反射率拟合结果。
23.具体地，一定大小的移动窗口m内，根据公式（3）利用建立基准期和预测期的降尺度数据间线性回归模型，利用最小二乘法求得移动窗口内中心像元的回归系数a和b，并应用于细影像，根据公式（4）估测得到的细影像；说明书其中，和分别为移动窗口m内和时期降尺度数据，为拟合残余误差，和分别为和时期细影像的反射率。
24.s3、残差的分配和补偿。采用平均聚合的方法将上述步骤2得到的的细影像的反射率升尺度至低空间分辨率，利用公式（5）计算其与粗影像反射率的差值，即粗分辨率下的残差rd；此项残差包含了上述步骤2中回归模型的拟合残差以及不同传感器间的误差；采用空间插值方法将残差分配至高空间分辨率，空间插值方法可以根据拟合情况选用适当的克里金插值法，将高空间分辨率下的残差rd补偿于的细影像，即细影像的反射率加上残差rd作为残差校正后的融合结果；式（5）中，表示低空间分辨率下的残差，和分别为
和时期粗像元的反射率，和分别为和时期粗像元中细像元的观测和估测反射率。
25.s4、拟合结果的空间滤波。在一定窗口内筛选中心像元的光谱相似像元，利用相似性像元的光谱差异性和与中心像元的相对说明书距离来加权计算预测期的中心像元值；具体地，根据公式（6），在w
×
w的窗口内，综合基准期细影像的各波段相似像元和中心像元的光谱差异d，差异最小的前s个像元被筛选为光谱相似像元；根据公式（7）计算窗口内邻近像元和中心像元的空间距离，并根据公式（8）确定各相似像元的权重；各光谱相似像元拟合结果的线性组合视为空间滤波结果。根据公式（9）通过加权窗口内相似性像元及其权重逐波段计算预测时期的中心像元值，实现对步骤3所得的估测反射率的空间滤波，从而得到时期空间滤波后的细影像；时期空间滤波后的细影像；时期空间滤波后的细影像；时期空间滤波后的细影像；其中，和分别表示细影像第b波段的i位置和中心位置像元的反射率，c为用于融合的影像波段数；w为窗口半径，为高空间分辨率下的残差分配值，为w
×
w窗口内中心像元的正向预测结果。
26.s5、正反向估测结果的加权：如果在时期之后可以获得第二说明书对粗-细影像，那么上述步骤1-4可以反向进行，正、反向估测结果的加权将有利于捕捉基准期与间发生的地表覆盖类别变化。设置起始日期（）和结束日期（）两个基准
期，由的粗-细影像估测的细影像称为正向估测，由的粗-细影像估测的细影像称为反向估测；采用线性加权方法综合正、反向估测结果作为加权时序融合结果；依据步骤s12得到的时期细影像的分类图，确定正、反向估测结果中每一地类的权重。具体地，以正向估测过程为例，由估测的细影像，根据公式（10），确定估测影像与实际观测影像间的回归系数k和l，并应用于的正向估测影像，根据公式（11）得到的似“真实影像”。根据公式（12）和（13）计算的正向估测影像和似“真实”影像的差距，得到正向估测的相对精度；向估测的相对精度；向估测的相对精度；向估测的相对精度；其中，和分别表示的正向估测和实际观测值，和分别表示的正向估测和似“真实”观测值，表示正向估测与实际观测值之间的差距，表示正向估测的相对精度；相似地，依照上述流程，由估测的细影像，确定估测影像说明书与实际观测影像间的回归系数，并应用于的反向估测影像，得到的似“真实影像”，计算的反向估测影像和似“真实”影像的差距，可以获得反向估测的相对精度；根据公式（14）和（15）确定各地类正、反向估测的权重；根据公式（16）获得正、反向估测加权的时序影像融合结果。向估测加权的时序影像融合结果。向估测加权的时序影像融合结果。
27.s6、结合高空间分辨率影像的连续校正：由于具有较高空间分辨率的sentinel-2、高分系列等卫星同时提供了较高的时间分辨率，即便受到天气条件等影响，通常也可以获
得多幅晴空下的影像，充分利用可获得的多幅细影像将有利于提升融合精度；基于上述步骤s5获得的加权时序影像融合结果，构建地表信息的动态变化模型；利用可获得的多幅已知细影像对动态模型进行修正，使动态模型模拟值与已知细影像值的误差最小，最终得到校正后的时空融合结果；具体地，采用s-g滤波方法对步骤s5得到的加权时序融合数据进行滤波去噪，并用于构建地表信息的动态模型；采用一种连续校正的cacao方法（公式（17）），利用多幅已知的细影像修正动态模型，通过最小化模型模拟值与实际观测值之间的误差，获得修正参数，将修正后的各波段反射率作为最终的时序遥感数说明书据融合结果；式（17）中，表示已知的细影像的各波段反射率，n为可获得的细影像的数量，表示动态模型模拟的各波段反射率，scale和shift为需要计算的修正参数，可通过最小二乘法获得；对于每一像元，可依据bppm*（ti+shift)获取任意时间的反射率值。
28.二、测试数据与处理本发明实施例以黑河流域中游大满超级站为例，所采用的卫星数据为sentinel-2多光谱影像（s2-msi）与准同步的中分辨率成像光谱仪影像所对应的6个波段，具体波段信息如表1所示。
29.modis数据选用2019年6月30日至2019年9月23日间逐天的mcd43a4产品，无云的sentinel-2影像时间分布如图2所示。
[0030] 说明书表1 s2-msiandmodis数据对应的波段信息
s2-msilevel-1影像从哨兵科学数据中心获取，作为细影像；说明书采用sen2cor(version2.8)进行大气校正，投影坐标为utmwgs-84，同时生成场景分类图，采用snap软件将各波段的空间分辨率重采样为20m；modis数据从美国航空局网站下载，作为粗影像；采用modis投影转换工具（mrt）重投影为utmwgs-84，空间分辨率重采样为460m。
[0031]
2019年6月30日作为基础期，2019年9月23日作为基础期，采用kmeans算法将基础期的sentinel-2影像分为10个类别；将此分类图用于modis影像的混合像元分解，窗口设置为31
×
31；利用基准期的modis降尺度影像和sentinel-2影像及预测期的modis降尺度影像分别进行正向和反向的估测，回归拟合中移动窗口内包含115
×
115个20m分辨率的像元，残差分配与空间滤波中窗口半径设为15，相似像元数为20；依据时期sentinel-2影像的分类图，获得每一地类在正、反向估测中的权重，并加权得到2019年6月30日-2019年9月23之间预测期的时间序列影像。此时间序列影像经过s-g滤波后构成地表信息动态模型，利用剩余日期的sentinel-2影像修正动态模型模拟结果，最终获得预测期连续修正后的高空间分辨率的时间序列遥感数据。
[0032]
三、结果验证采用留一影像交叉验证的方法评估融合影像的质量，即在连续校正过程中，预留一期的sentinel-2影像用于结果验证，其余日期的影像用于修正动态模型，依次对比各日期的融合影像与实说明书际观测影像；定性评价采用空间可视化图评估融合影像的空间特征；定量评价选用融合影像与实际观测影像间的相关系数cc、质量指数qi、均方根误差rmse以及线性拟合系数slope和截距intercept，计算公式如下：
(18)(19)(20)cc用于评价估测影像反射率与实际观测影像反射率之间的相关性，负数表示负相关，正数表示正相关。rmse表示融合误差，值越小表示误差越小；qi用于评价影像的结构相似性，值越接近1表示越相似。
[0033]
下面以2019年8月4日为例，对比原fit-fc方法与本发明方法在研究区内对sentinel-2和modis影像时空融合的表现。表2为原fit-fc方法与本发明方法估测影像的各波段反射率与实际观测反射率的定量评价指标，图3为两种方法融合影像和实际观测影像各波段反射率的散点密度图，图4为两种方法融合影像与实际观测影像的空间可视化图；可以看出，本发明方法融合生成的影像与实际观测影像更接近，消除了原算法产生的模糊现象，呈现出更丰富的纹理特征，极大地提高了融合精度指标，可以提供更准确的反射率。
[0034]
说明书表2各波段反射率与实际观测反射率的定量评价指标