一种高寒地区长时间序列高频率NDVI生成方法和装置

一种高寒地区长时间序列高频率ndvi生成方法和装置
技术领域
1.本发明属于遥感技术领域，尤其涉及一种高寒地区长时间序列高频率ndvi生成方法和装置。


背景技术：

2.归一化植被指数(ndvi)是植被生长状态及植被覆盖度的最佳指示因子。许多研究表明植被生长周期内的ndvi对于降水量、co2浓度的季节性变化和地理位置的变化均有反映。因此，ndvi可以作为监测植被和生态环境变化的有效指标。
3.现有长时间序列遥感影像数据主要有modis影像和landsat影像，modis影像时间分辨率高，在植被生长期内可用影像较多，但空间分辨率不足，不适用于小范围区域内ndvi计算；landsat有较高的空间分辨率，可以满足小范围区域内ndvi的计算要求，但其时间分辨率较低，每月仅有2景影像。但高寒地区只有7
‑
9月较为短暂的植被生长期，且受多云雨天气影响，常出现可用landsat影像受限、高质量影像缺失的问题，导致无法计算长时间序列高频率ndvi，无法对植被覆盖情况进行高时间频率监测。为实现高寒地区长时间序列高频率ndvi的计算，亟需一种适用于高寒地区长时间序列高频率ndvi的计算方法。


技术实现要素：

4.本发明要解决的技术问题是，提供一种高寒地区长时间序列高频率ndvi生成方法和装置，以克服针对高寒地区植被生长期短、可使用影像时间受限、高质量遥感影像缺失，无法计算长时间序列ndvi的缺陷。
5.为实现上述目的，本发明采用如下的技术方案：
6.一种高寒地区长时间序列高频ndvi生成方法，包括以下步骤：
7.步骤1、对于t2时刻缺失landsat
‑
8影像的情况，获取t1和t2时刻低空间分辨率高时间分辨率的modis数据和t1时刻高空间分辨率低时间分辨率的landsat
‑
8数据；
8.步骤2、将步骤1中的t1和t2时刻低空间分辨率高时间分辨率的modis数据以及t1时刻高空间分辨率低时间分辨率的landsat
‑
8数据利用亚像元类分数变化信息的融合方法进行时空融合，得到t2时刻的landsat
‑
8数据；
9.步骤3、对于t2时刻缺失landsat
‑
8影像的情况，采用sentinel
‑
2卫星影像填补缺失的landsat
‑
8影像；
10.步骤4、对于t2时刻存在landsat
‑
8影像但有云阴影覆盖的情况，利用邻域相似像素插值方法将含云的landsat
‑
8影像进行去云处理，得到去云后的结果；
11.步骤5、对步骤2中时空融合的结果、步骤3中sentinel
‑
2卫星影像和步骤4影像去云的结果分别进行影像矫正；
12.步骤6、根据所述三个影像矫正结果分别计算ndvi。
13.作为优选，步骤2具体为：
14.步骤21、估计t1时刻精细图像的端元与类丰度；
15.步骤22、估计t2时刻精细图像的端元与类别分数；
16.步骤23、根据类丰度和类别分数确定从t1到t2的端元变化和亚像元覆盖率变化信息，同时对t2处的精细图像进行时间预测；
17.步骤24、采用空间预测策略来修正所述时间预测，得到sfsdaf的预测结果，即t2时刻的landsat
‑
8数据。
18.作为优选，步骤4具体为：
19.步骤41、采用自适应的移动窗口搜索方法来选取邻域相似像素；
20.步骤42、根据相似像素的权重，并且利用目标图像中所有相似像素的加权平均，对目标像素进行预测；
21.步骤43、根据预测结果，将含云的landsat
‑
8影像进行去云处理，得到去云后的结果。
22.作为优选，对时空融合结果进行矫正具体包括：
23.首先将t1时刻时空融合的结果与t1时刻原始影像利用fcm算法做变化检测，得到变化检测结果，将变化区域和未变化区域分别与原始影像建立线性模型；然后再将t2时刻时空融合的结果与t2时刻原始影像利用fcm算法做变化检测，得到变化检测结果；最后将在t1时刻求得的线性模型分别应用于t2时刻时空融合结果的变化区域和未变化区域，预测t2时刻的landsat
‑
8影像，实现对t2时刻时空融合结果的矫正。
24.作为优选，对sentinel
‑
2卫星影像进行矫正具体为：
25.首先将t1时刻sentinel
‑
2卫星影像与t1时刻原始影像建立线性模型，然后再将在t1时刻求得的线性模型应用于t2时刻sentinel
‑
2卫星影像，预测t2时刻的landsat
‑
8影像，实现对t2时刻sentinel
‑
2卫星影像的矫正。
26.作为优选，对含云影像去云处理后的结果进行矫正具体为：
27.将无云区域的像元与原始影像对应区域建立线性模型，求得模型后将其应用于有云覆盖并去云后的区域，实现对影像去云结果的矫正。
28.本发明还提供一种高寒地区长时间序列高频ndvi生成装置，包括：
29.第一处理模块，用于对于t2时刻缺失landsat
‑
8影像的情况，对t1和t2时刻低空间分辨率高时间分辨率的modis数据和t1时刻高空间分辨率低时间分辨率的landsat
‑
8数据，利用亚像元类分数变化信息的融合方法进行时空融合，得到t2时刻的landsat
‑
8数据；
30.第二处理模块，用于对于t2时刻缺失landsat
‑
8影像的情况，采用sentinel
‑
2卫星影像填补缺失的landsat
‑
8影像；
31.第三处理模块，用于对于t2时刻存在landsat
‑
8影像但有云阴影覆盖的情况，利用邻域相似像素插值方法将含云的landsat
‑
8影像进行去云处理，得到去云后的结果；
32.矫正模块，用于对时空融合的结果、sentinel
‑
2卫星影像和影像去云的结果分别进行影像矫正；
33.计算模块，用于根据所述三个影像矫正结果分别计算ndvi。
34.作为优选，第一处理模块包括：
35.第一估计单元，用于估计t1时刻精细图像的端元与类丰度；
36.第二估计单元，用于估计t2时刻精细图像的端元与类别分数；
37.第一预测单元，用于根据类丰度和类别分数确定从t1到t2的端元变化和亚像元覆
盖率变化信息，同时对t2处的精细图像进行时间预测；
38.修正单元，用于采用空间预测策略来修正所述时间预测，得到sfsdaf的预测结果，即t2时刻的landsat
‑
8数据。
39.作为优选，第三处理模块包括：
40.搜索单元，用于采用自适应的移动窗口搜索方法来选取邻域相似像素；
41.第二预测单元，用于根据相似像素的权重，并且利用目标图像中所有相似像素的加权平均，对目标像素进行预测；
42.去云单元，用于根据预测结果，将含云的landsat
‑
8影像进行去云处理，得到去云后的结果。
43.本发明的高寒地区长时间序列高频率ndvi生成方法和装置，利用时空融合和影像去云相结合方式来弥补缺失遥感影像，从而实现长时间序列高频率ndvi的计算。
附图说明
44.图1为本发明高寒地区长时间序列高频率ndvi生成方法流程图；
45.图2为本发明高寒地区长时间序列高频率ndvi生成装置的结构示意图。
具体实施方式
46.如图1所示，本发明提供一种高寒地区长时间序列高频率ndvi生成方法，其特征在于包括以下步骤：
47.步骤1：数据获取：针对t2时刻缺失landsat
‑
8影像的情况，获取t1和t2时刻低空间分辨率且高时间分辨率的modis数据和t1时刻高空间分辨率且低时间分辨率的landsat
‑
8数据。
48.步骤2：时空融合：将步骤1中的三幅影像利用亚像元类分数变化信息的融合方法(sfsdaf)进行时空融合，得到t2时刻的landsat
‑
8数据。
49.步骤3：针对t2时刻缺失landsat
‑
8影像的情况，采用sentinel
‑
2卫星影像填补缺失的landsat
‑
8影像。
50.步骤4：影像去云：针对t2时刻存在landsat
‑
8影像但有云阴影覆盖的情况，利用邻域相似像素插值(neighborhood similar pixel interpolation,nspi)将含云的landsat
‑
8影像进行去云处理，得到去云后的结果。
51.步骤5：结果矫正：对步骤2中时空融合的结果、步骤3中sentinel
‑
2卫星影像和步骤4影像去云的结果分别使用不同的方法进行线性拟合，实现影像矫正并计算其与原始影像间的均方根误差和相关系数。
52.步骤6：计算ndvi：利用步骤5中矫正的结果，分别计算ndvi。
53.进一步，步骤2中采用sfsdaf的方法进行时空融合，该方法不仅直接推导端元变化来表示物候变化，而且还直接在精细像素中推导亚像元土地覆盖类分数变化以适应土地覆盖类型的变化以及精细图像比例下是否存在混合像元。该方法有四个步骤：
54.步骤21、估计t1时刻精细图像的端元与类丰度。先使用非监督分类算法(k
‑
means)将t1时刻landsat
‑
8影像fr聚类为土地覆盖图，在生成土地覆盖图之后，每个端元的反射率是聚类图中像素的平均值。非监督分类产生的是一种硬分类图，在sfsdaf中，需要生成软分
类图以指示亚像元尺度的土地覆盖信息，类丰度计算如下：
[0055][0056]
其中，v(x
ij
，y
ij
，t1)是t1时刻fr中(x
ij
，y
ij
)处的反射向量，μ(c，t1)是t1时刻第c个簇中心。
[0057]
步骤22、估计t2时刻精细图像端元与类别分数。首先对计算了两个不同时间点粗略图像的类别分数，并得到了在研究时间段内地表事物在粗略图像尺度上的类别分数变化。然后利用双三次插值法将粗略尺度的类别分数变化插值到精细尺度上，并使用了基于空间距离的加权策略来修正所得到的精细尺度的类别分数变化。之后利用粗略图像的时间差异和类别分数来计算端元的变化。
[0058]
步骤23、同时考虑从t1到t2的端元变化和亚像元覆盖率变化信息，对t2处的精细图像进行时间预测。
[0059][0060]
其中，表示t1处第b波段的fr反射率，表示t2时刻fr图像地物覆盖率与端元的反射率，则对应的表示t1时刻。
[0061]
步骤24、与经典的融合算法fsdaf相似，引入空间预测策略来修正步骤23产生的时间预测，就得到了sfsdaf的最终预测结果。
[0062][0063]
其中，表示t1时刻第b波段的fr反射率，δr
fr
(x
k
，y
k
，b)表示fr总反射率变化。
[0064]
进一步，步骤4中利用nspi方法对含云影像进行去云处理时，采用自适应的移动窗口搜索方法来选取邻域相似像素。将含云的影像定义为目标影像，而选择填充目标图像中含云区域的其他图像称为输入影像。
[0065]
首先，获取与目标影像日期接近的landsat
‑
8影像作为输入影像；
[0066]
然后，采用自适应的移动窗口搜索方法来选择邻近相似像元。假设在含云区域附近的相同土地覆盖像元与目标缺失像元具有相似的光谱特征和时间变化规律，根据光谱相似性从共同像素中选择相似像素。将光谱相似性定义为每个公共像素(目标和输入图像中均为有效值的含云区域外的所有公共像素)与目标像素(位于目标图像含云区域中没有有效值的像素)之间的均方根偏差(rmsd)。
[0067]
[0068]
其中，l(x
i
,y
i
,t1,b)是t1时刻获取的输入图像在b波段(x
i
,y
i
)处的第i个公共像素的值，l(x,y,t1,b)是目标图像的对应值，n是光谱波段数。
[0069]
然后，使用一个阈值来识别rmsd值低于阈值的类似像素。阈值可以由输入图像的像素群的标准差和估计的土地覆盖类别的数量来确定。如果第i个公共像素的rmsd满足下式，则选择第i个公共像素作为相似像素：
[0070][0071]
其中，σ(b)为波段的整个输入图像的标准差，m为类数。估计的类数量(m)需要预定义。这个值是一个经验阈值，随着景观的复杂性而变化。它可以通过对输入图像的目视解译来估计，或者使用先前的土地覆盖地图。在这个实验中，本发明使用m的值为5。
[0072]
然后，计算相似像素的权重。所有相似像素点的信息可以用来预测目标像素点的值。然而，相似像素的贡献可能会有所不同，因为一些相似像素比其他相似像素更有可能与目标像素在光谱上具有可比性。权重w
j
决定了第j个相似像素对预测目标像素值的贡献。这是由相似像素的位置和相似像素与目标像素之间的光谱相似性决定的。光谱相似性越高，与目标像素的距离越小，该像素的权值越高。计算第j个相似像素(x
j
,y
j
)与目标像素(x,y)之间的地理距离d
j
：
[0073][0074]
光谱相似度由每个相似像元与目标像元之间的rmsd决定。结合光谱相似性和地理距离，综合指数cd可计算为：
[0075]
cd
j
＝rmsd
j
×
d
j
[0076]
相似像素cd值越大，对目标像素的计算值的贡献就越小，所以用cd的归一化倒数作为权重w
j
:
[0077][0078]
最后，计算目标像素值。每个相似像素所提供的信息的可靠性可能是不同的，相似像素的权重越大，越可靠。据此，利用目标图像中所有相似像素的加权平均，对目标像素进行预测：
[0079][0080]
进一步，步骤5结果矫正中对步骤2中时空融合的结果、步骤3中sentinel
‑
2卫星影像和步骤4影像去云的结果分别进行矫正。
[0081]
(1)对时空融合结果进行矫正。将t1时刻时空融合的结果与t1时刻原始影像利用fcm算法做变化检测，得到变化检测结果，将变化区域和未变化区域分别与原始影像建立线性模型。然后再将t2时刻时空融合的结果与t2时刻原始影像利用fcm算法做变化检测，得到变化检测结果，最后将在t1时刻求得的线性模型分别应用于t2时刻时空融合结果的变化区域和未变化区域，从而预测t2时刻的landsat
‑
8影像，实现对t2时刻时空融合结果的矫正，将
矫正结果与t2时刻原始影像作对比，计算其均方根误差和相关系数。
[0082]
(2)对sentinel
‑
2卫星影像进行矫正。利用sentinel
‑
2卫星影像填补缺失的landsat
‑
8影像，由于是异源遥感影像，要对其进行传感器间的矫正。将t1时刻sentinel
‑
2卫星影像与t1时刻原始影像建立线性模型，然后再将在t1时刻求得的线性模型应用于t2时刻sentinel
‑
2卫星影像，从而预测t2时刻的landsat
‑
8影像，实现对t2时刻sentinel
‑
2卫星影像的矫正，将矫正结果与t2时刻原始影像作对比，计算其均方根误差和相关系数。
[0083]
(3)对含云影像去云处理后的结果进行矫正。将无云区域的像元与原始影像对应区域建立线性模型，求得模型后将其应用于有云覆盖并去云后的区域，实现对影像去云结果的矫正，最后计算影像去云并矫正的结果与原始影像间的均方根误差和相关系数。
[0084]
如图2所示，本发明还提供一种高寒地区长时间序列高频ndvi生成装置，包括：
[0085]
第一处理模块，用于对于t2时刻缺失landsat
‑
8影像的情况，对t1和t2时刻低空间分辨率高时间分辨率的modis数据和t1时刻高空间分辨率低时间分辨率的landsat
‑
8数据，利用亚像元类分数变化信息的融合方法进行时空融合，得到t2时刻的landsat
‑
8数据；
[0086]
第二处理模块，用于对于t2时刻缺失landsat
‑
8影像的情况，采用sentinel
‑
2卫星影像填补缺失的landsat
‑
8影像；
[0087]
第三处理模块，用于对于t2时刻存在landsat
‑
8影像但有云阴影覆盖的情况，利用邻域相似像素插值方法将含云的landsat
‑
8影像进行去云处理，得到去云后的结果；
[0088]
矫正模块，用于对时空融合的结果、sentinel
‑
2卫星影像和影像去云的结果分别进行影像矫正；
[0089]
计算模块，用于根据所述三个影像矫正结果分别计算ndvi。
[0090]
进一步，第一处理模块包括：
[0091]
第一估计单元，用于估计t1时刻精细图像的端元与类丰度；
[0092]
第二估计单元，用于估计t2时刻精细图像的端元与类别分数；
[0093]
第一预测单元，用于根据类丰度和类别分数确定从t1到t2的端元变化和亚像元覆盖率变化信息，同时对t2处的精细图像进行时间预测；
[0094]
修正单元，用于采用空间预测策略来修正所述时间预测，得到sfsdaf的预测结果，即t2时刻的landsat
‑
8数据。
[0095]
进一步，第三处理模块包括：
[0096]
搜索单元，用于采用自适应的移动窗口搜索方法来选取邻域相似像素；
[0097]
第二预测单元，用于根据相似像素的权重，并且利用目标图像中所有相似像素的加权平均，对目标像素进行预测；
[0098]
去云单元，用于根据预测结果，将含云的landsat
‑
8影像进行去云处理，得到去云后的结果。
[0099]
实施例1：
[0100]
采用2020年不同时期获取的modis和landsat
‑
8卫星数据进行时空融合的实验，sentinel
‑
2卫星数据用来填补缺失的landsat
‑
8影像，数据详细信息如表1所示，采用遥感影像的红波段和近红外波段进行ndvi的计算，研究区域为扎赉诺尔灵泉露天煤矿。
[0101]
表1遥感影像数据信息
[0102][0103]
本发明为高寒地区长时间序列高频ndvi生成方法——以扎赉诺尔灵泉露天煤矿为例，包括以下步骤：
[0104]
步骤1：数据获取。针对t2时刻缺失landsat
‑
8影像的情况，获取t1和t2时刻低空间分辨率高时间分辨率的modis数据和t1时刻高空间分辨率低时间分辨率的landsat
‑
8数据。
[0105]
步骤2：时空融合。将步骤1中的三幅影像利用亚像元类分数变化信息的融合方法(sfsdaf)进行时空融合，得到t2时刻的landsat
‑
8数据。该方法首先在t1时刻估计精细图像的端元与类丰度；然后估计t2处的精细图像端元与类别分数；同时考虑从t1到t2的端元变化和亚像元覆盖率变化信息，对t2处的精细图像进行时间预测；最后引入空间预测策略来修正上述预测，从而得到sfsdaf的最终预测结果。
[0106]
步骤3：针对t2时刻缺失landsat
‑
8影像的情况，采用sentinel
‑
2卫星影像填补缺失的landsat
‑
8影像。
[0107]
步骤4：影像去云。针对t2时刻存在landsat
‑
8影像但有云阴影覆盖的情况，利用邻域相似像素插值(neighborhood similar pixel interpolation,nspi)，采用自适应的移动窗口搜索方法来选取邻域相似像素，通过计算相似像素的权重，利用目标图像中所有相似像素的加权平均，对目标像素进行预测，最终将含云的landsat
‑
8影像进行去云处理，得到去云后的结果。
[0108]
步骤5：结果矫正。对步骤2中时空融合的结果、步骤3中sentinel
‑
2卫星影像和步骤4影像去云的结果分别使用不同的方法与原始影像进行线性拟合。
[0109]
最后，计算其与原始影像间的均方根误差和相关系数，结果如表2所示，根据相关系数的评判标准，大于0.8认为高度相关，因此，可以一定程度上解决可用landsat影像受限，高质量影像缺失，导致无法计算长时间序列高频率ndvi的问题。
[0110]
表2矫正结果分析
[0111] 时空融合sentinel
‑
2影像去云均方根误差(rmse)0.13280.10740.1153相关系数(r)0.88640.93570.8996
[0112]
步骤6：计算ndvi。利用步骤5中矫正的结果，分别计算ndvi。
[0113]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。