一种基于时空加权的生产高时空分辨率NDVI的方法与流程

本发明涉及一种多源遥感数据融合的方法，主要用来生产高时空分辨率的归一化植被指数(ndvi)，属于时空数据融合的范畴。

背景技术：

归一化植被指数(normalizeddifferencevegetationindex，缩写为ndvi)是一种广泛使用的植被指数。ndvi数值能够反映植被的生长状况并与植被的生理参数比如叶绿素含量，含水量以及叶面积指数等有直接的联系。

cn102831310b公开了一种构建高空间分辨率的时间序列数据的方法，其背景技术部分对于ndvi数据的相关知识进行了详细说明，其中提及两种现有的ndvi数据：landsatndvi数据和modisndvi数据。其中，landsatndvi数据中的ndvi数值的空间分辨率为30m*30m，时间分辨率为16天，而modisndvi数据中的ndvi数值的最大空间分辨率为250m*250m，时间分辨率为1天。也就是说，landsatndvi数据具有较高的空间分辨率和较低的时间分辨率，而modisndvi数据具有较低的空间分辨率和较高的时间分辨率。

为了能够同时获得高空间分辨率和高时间分辨率的ndvi数据，现有技术发展出了多种对landsatndvi数据和modisndvi数据进行融合，以获取高空间分辨率的ndvi数据的方法，主要包括三种类型：(1)基于权函数的方法；(2)基于混合像元分解的方法；(3)基于混合模型的方法。

其中，(1)基于权函数的方法以starfm(空间和时间自适应反射融合模型，spatialandtemporaladaptivereflectancefusionmodel)为典型代表，参见gao，f.，masek，j.，schwaller，m.，&hall，f.(2006)，“对landsat与modis的地表反射率进行融合：预测每日的landsat地表反射率”，电气与电子工程师学会之地球科学与遥感学报，44，2207-2218。(gao，f.，masek，j.，schwaller，m.，&hall，f.(2006).ontheblendingofthelandsatandmodissurfacereflectance:predictingdailylandsatsurfacereflectance.ieeetransactionsongeoscienceandremotesensing，44，2207-2218)。该类方法着眼于相似像元的搜索及其权重的确定，通过当前像元与周围像元在时间、空间以及光谱等特征上的相似性确定相似像元并计算对应的权重，以相似像元的平均增量作为当前像元的增量。这类方法往往在纯净度比较高的区域才有比较好的效果。

(2)基于混合像元分解的方法以ndvi-lmgm(ndvi-线性混合生长模型，ndvilinearmixinggrowthmodel)为代表，参见rao，y.，zhu，x.，chen，j.，&wang，j.(2015)，“一种升级过地使用多时相modisndvi数据与landsattm/etm+影像生产高时空分辨率的ndvi时序数据的方法”，遥感，7，7865-7891(rao，y.，zhu，x.，chen，j.，&wang，j.(2015).animprovedmethodforproducinghighspatial-resolutionndvitimeseriesdatasetswithmulti-temporalmodisndvidataandlandsattm/etm+images.remotesensing，7，7865-7891)。该类方法通过混合像元分解，计算出每个像元在时间上的增量，从而得到预测时间的ndvi数据，有关该方法的具体内容也可参见前述现有技术cn102831310b，二者实质上是相同的。这类方法对于异质性程度比较高的区域也能够有比较好的效果，但其假设土地覆盖状况不发生变化，这一点不一定能够满足。

(3)基于混合模型的方法以fsdaf(灵活时空数据融合方法，flexiblespatiotemporaldatafusionmethod)为代表，参见zhu，x.，helmer，e.h.，gao，f.，liu，d.，chen，j.，&lefsky，m.a.(2016)，“一种灵活地将不同分辨率的卫星影像进行时空融合的方法”，环境遥感，172，165-177(zhu，x.，helmer，e.h.，gao，f.，liu，d.，chen，j.，&lefsky，m.a.(2016).aflexiblespatiotemporalmethodforfusingsatelliteimageswithdifferentresolutions.remotesensingofenvironment，172，165-177)。该类方法将多种模型，如权函数模型、混合像元分解模型以及空间插值模型等结合起来，计算得到高分辨率影像上每个像元的增量。这类方法往往比较复杂，但是精度相对较高，而且能够对土地覆盖变化进行一定程度的捕捉。

以上三类方法在各自适用的领域里都取得了比较好的效果。但是它们都没有考虑到使用两期低空间分辨率数据插值到高空间分辨率的方式，两期插值结果之间的增量能够从一定程度上反映像元尺度的增量，同时也包含了土地覆盖变化的信息。本发明基于贝叶斯模型平均方法将混合像元分解的增量与空间插值的增量进行综合，能够有效提高数据融合的精度。为生产高时空分辨率的ndvi提供了新的途径。

技术实现要素：

本发明要解决的技术问题是提供一种基于时空加权的生产高时空分辨率ndvi的方法，以减少或避免前面所提到的问题。

为解决上述技术问题，本发明提出了一种基于时空加权的生产高时空分辨率ndvi的方法，用于将高空间分辨率和低时间分辨率的landsatndvi数据与低空间分辨率和高时间分辨率的modisndvi数据融合，从基准时间t1获得预测时间tp的高时空分辨率的ndvi数据，所述方法包括如下步骤：

步骤a：以同一基准时间t1为起点，利用modisndvi数据和landsatndvi数据，分别获得预测时间tp的每个像元在ndvi上的第一增量△ndvid和第二增量△ndvis；其中，

所述第一增量△ndvid指的是，从基准时间t1到预测时间tp，假定modisndvi数据的每个像元的ndvi增量，等于同样面积下的landsatndvi数据的不同地物类别的像元的ndvi增量的平均值，则利用t1时间的modisndvi数据和landsatndvi数据，以及tp时间的modisndvi数据，计算获得tp时间到t1时间的landsatndvi数据的每个像元的ndvi增量，该增量即为所述第一增量△ndvid；

所述第二增量△ndvis指的是，采用现有envi/idl软件中集成的薄板样条方法工具，直接将t1时间与tp时间的低空间分辨率modis像元分别插值到高空间分辨率的landsat像元上去，然后使用tp时间的插值结果与t1时间的插值结果相减，其结果即为每个像元的所述第二增量△ndvis；

步骤b，使用贝叶斯模型平均方法，将步骤a获得的每个高空间分辨率的像元在ndvi上的所述第一增量与所述第二增量进行权重组合计算，获得每个像元的ndvi综合增量；

步骤c，通过步骤b计算得到的t1时间到tp时间，每个像元的所述ndvi综合增量，加到t1时间对应的高空间分辨率ndvi数值上，即可得到tp时间的高空间分辨率ndvi数值。

优选地，所述步骤b中，每个高空间分辨率像元(xi，yj)的所述ndvi综合增量用公式表示为，

δndvicom(xi,yj)＝ws*δndvis(xi,yj)+wd*δndvid(xi,yj)(2)

其中，公式(2)中，wd和ws分别表示所述第一增量和所述第二增量的权重，△ndvicom(xi，yj)表示所述第一增量与所述第二增量结合得到的所述ndvi综合增量。

优选地，所述公式(2)种的每个像元的所述权重wd和ws通过贝叶斯模型平均方法求解获得。

优选地，所述步骤b进一步包括对所述ndvi综合增量进行误差修正的步骤。

优选地，所述误差修正的步骤包括：

首先计算出低空间分辨率像元(x，y)整体的误差为，

公式(4)中，m表示低空间分辨率像元(x，y)所覆盖的高空间分辨率像元的数目；r(x，y)表示低空间分辨率像元(x，y)整体的误差；△ndvi^c(x，y)表示tp时间与t1时间低空间分辨率像元(x，y)的ndvi的差值；

然后计算高空间分辨率像元的均质性如下，

hi(xi,yj)＝nsmae/m，(5)

公式(5)中，hi(xi，yj)是以高空间分辨率像元(xi，yj)为中心，构建的8*8的窗口的均质性参数；m表示这个窗口中高空间分辨像元的数目，nsame表示这个窗口中，与中心像元(xi，yj)具有相同类型的高空间分辨率像元的数目；

之后，根据像元的均质性来对误差进行分配如下，

公式(6)中，r(xi，yj)即为高空间分辨率像元(xi，yj)应该分配到的误差；

最后，将这个误差加到公式(2)中的综合增量上去，得到修正过的高空间分辨率像元的ndvi增量如下，

公式(7)中，△ndvicom^new(xi，yj)即为最后计算得到的高空间分辨率像元(xi，yj)从t1时间到tp时间的修正过的所述ndvi综合增量。

优选地，所述步骤c中，利用所述修正过的ndvi综合增量计算tp时间的高空间分辨率ndvi的公式为：

公式(8)中，ndvip(xi，yj)即为tp时间的每个像元(xi，yj)的高空间分辨率ndvi，ndvi1(xi，yj)即为t1时间的每个像元(xi，yj)的高空间分辨率ndvi。

本发明通过对低空间分辨率的数据和高空间分辨率数据进行混合像元分解以及空间插值得到高空间分辨率数据在时间上的两种增量，然后采用贝叶斯模型平均方法将所述两种增量进行权重组合计算，得到优化过的综合增量，提高了增量计算的准确性，充分利用了地理数据间的时空联系，提高了融合精度。

附图说明

以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范围。其中，

图1显示的是根据本发明的一个具体实施例的第一增量的原理示意图；

图2a-2f分别显示的是土地覆盖变化区域以及对应的landsatndvi数据和modisndvi数据；其中表示的某集水流域区的情况，分别为：图2a表示2004年11月26日的landsatndvi，图2b表示的2004年12月12日landsatndvi；从上述图像分别模拟得到图2d表示的2004年11月26日和图2e表示的2004年12月12日modisndvi；图2c表示的是2004年11月26日的landsat影像，图2f表示对应的非监督分类结果。

图3a-3e分别显示的是几种不同的数据融合方法在土地覆盖变化区域的融合结果；其中图3a-3d分别表示的是2004年12月12日landsatndvi经过ndvi-lmgm融合、starfm融合、fsdaf融合、本发明的方法融合的结果，图3e表示的是真实值。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图说明本发明的具体实施方式。但本领域的技术人员应该知道，以下实施例并不是对本发明技术方案作的唯一限定，凡是在本发明技术方案精神实质下所做的任何等同变换或改动，均应视为属于本发明的保护范围。

下面详细说明根据本发明提供的基于时空加权的生产高时空分辨率ndvi的方法的原理，当然，与背景部分所述的三种类型的对landsatndvi数据和modisndvi数据进行融合的方法的一样，本发明也是对上述不同数据源的遥感数据进行融合的方法，亦即，本发明的基于时空加权的生产高时空分辨率ndvi的方法用于将高空间分辨率和低时间分辨率的landsatndvi数据与低空间分辨率和高时间分辨率的modisndvi数据融合，从基准时间t1获得预测时间tp的高时空分辨率的ndvi数据。

由于获取landsatndvi数据和modisndvi数据的卫星有相似的轨道参数以及不足30分钟的卫星过境时间间隔(该间隔对于每天24小时来说可以忽略不计)，因此，在理论上可以假设这两种数据同一地点和同一基准时间t1的ndvi数据仅存在空间分辨率上的差别，其余均可忽略不计。其中，landsatndvi数据的空间分辨率为30m*30m，时间分辨率为16天；modisndvi数据的空间分辨率为250m*250m，时间分辨率为1天；则数据融合之后生产获得的高时空分辨率的ndvi数据，其最大空间分辨率可以等于landsatndvi数据的空间分辨率(30m*30m)，其最大时间分辨率可以等于modisndvi数据的时间分辨率(1天)。也就是说，通过已知空间分辨率为250m*250m、时间分辨率为1天的modisndvi数据，以及空间分辨率为30m*30m、时间分辨率为16天的landsatndvi数据，以同一基准时间t1为起点，本发明的方法可以预测获得其后15天的任意1天的预测时间tp的高时空分辨率的ndvi数据。

下面详细说明本发明的基于时空加权的生产高时空分辨率ndvi的方法的具体步骤。

步骤a：以同一基准时间t1为起点，利用modisndvi数据和landsatndvi数据，分别获得预测时间tp的每个像元在ndvi上的第一增量△ndvid和第二增量△ndvis。

其中，所述第一增量△ndvid指的是，从基准时间t1到预测时间tp，假定modisndvi数据的每个像元的ndvi增量，等于同样面积下的landsatndvi数据的不同地物类别的像元的ndvi增量的平均值，则利用t1时间的modisndvi数据和landsatndvi数据，以及tp时间的modisndvi数据，计算获得tp时间到t1时间的landsatndvi数据的每个像元的ndvi增量，该增量即为所述第一增量△ndvid。

具体来说，基准时间t1(实际上存在30分钟的时间间隔)的低空间分辨率modisndvi数据和高空间分辨率的landsatndvi数据对同一地点来说仅存在空间分辨率上的差别，为了融合时将两种数据的位置对准方便计算，可以将低空间分辨率的modisndvi数据从250m*250m重采样到240*240m，从而一个modis像元可以完全覆盖8*8＝64个landsat像元。

如图1所示，其显示的是根据本发明的一个具体实施例的第一增量的原理示意图，其中时间轴的上方表示的是低空间分辨率的modisndvi数据，时间轴下方表示的是高空间分辨率的landsatndvi数据。以同一基准时间t1为起点，利用t1时间的modisndvi数据和landsatndvi数据，以及tp时间的modisndvi数据，可以获得时间轴右下方虚线表示的原本不存在的tp时间的高空间分辨率的landsatndvi数据。

此时，假定从基准时间t1到预测时间tp，modisndvi数据的每个像元的ndvi增量，等于同样面积下的landsatndvi数据的不同地物类别的像元的ndvi增量的平均值。即，图1中时间轴上方心形表示的1个modis像元从t1至tp时间的ndvi增量，可以通过将tp时间的低空间分辨率的modisndvi数据与t1时间的低空间分辨率的modisndvi数据直接相减得到。

对应于时间轴上方的1个modis像元的面积，等于时间轴下方的64个landsat像元的面积，将这64个landsat像元的ndvi增量按照其面积平均起来，即可认定为等于时间轴上方的同样面积的1个modis像元的ndvi增量。其计算过程为，根据高空间分辨率的landsatndvi数据的地物分类计算出每个低空间分辨率的modisndvi数据的像元中每种地物类型的丰度。分类数据可以来源于已有的分类产品，比如globeland30，也可以通过对高空间分辨率的landsatndvi数据的卫星影像数据进行分类得到。根据高空间分辨率的landsatndvi数据的地物分类结果，计算出每个低空间分辨率的modisndvi数据的像元中包含的不同地物类型的面积比例，即丰度。因此，将每种地物类型的ndvi增量与其丰度的乘积相加，就可以得到该面积下的全部landsat像元(64个像元)的ndvi增量的平均值，这个平均值正好等于前述同样面积的modis像元(1个像元)从t1至tp时间的ndvi增量。

例如，图1时间轴下方的64个landsat像元包括三种地物类型，则第1类地物类型的24个像元具有一个ndvi增量，第2类地物类型的20个像元具有另一个ndvi增量，第3类地物类型的20个像元具有又一个ndvi增量，这三种ndvi增量，相对于每个像元来说都是其从t1至tp时间的第一增量△ndvid。为了计算出这三个增量，需要构建至少三个方程，此时假定最邻近时间轴上方的心形标记像元的8个像元(亦即围绕心形标记像元的8个像元)的ndvi增量与该心形标记像元的ndvi增量相同，则可以以每个低空间分辨率像元为中心，构建3*3的局部窗口。假设这个局部窗口中每一种地物类型在t1时间至tp时间内的ndvi增量是相同的。根据这个局部窗口中的9个低空间分辨率像元的ndvi增量，可构建如下的线性方程组，公式表示如下，

公式(1)中，△ndvi^c(x，y)表示从t1到tp时间低空间分辨率像元(x，y)的ndvi增量，fl(x，y)表示第l种类型在低空间分辨率像元(x，y)中的丰度，△ndvil(x，y)表示3*3的局部窗口中第l种类型的ndvi增量。

根据公式(1)通过最小二乘进行求解可以得到△ndvil(x，y)。对于图示局部窗口中的低空间分辨率像元(x，y)而言，其覆盖范围下的高分辨率像元(xi，yj)在时间上的增量为△ndvid(xi，yj)，这个增量即称为第一增量，下标(i，j)表示低空间分辨率像元覆盖下的第(i，j)个高空间分辨率像元。已经假设这个局部窗口中，每种地物类型的ndvi增量是固定的。因此，对于高空间分辨率像元(xi，yj)而言，其在时间上的增量△ndvid(xj，yj)就等于高空间分辨率像元(xj，yj)所属的那个地物类型的ndvi增量。如果高空间分辨率像元(xi，yj)属于第l类，那么△ndvid(xi，yj)就等于△ndvil(x，y)。

对于第一增量△ndvid而言，由于tp时间的地物类型与t1时间的地物类型可能有较大不同，例如水土流失、森林砍伐、火灾、洪水等均可造成不同时间的地物类型发生变化。因此，公式(1)中的每种地物类型的丰度fl(x，y)在两个时间下有可能是不同的，会导致地物类型的ndvi增量△ndvil(x，y)计算出错。而且在将每种地物类型的ndvi增量△ndvil(x，y)按类型分配给tp时间的高空间分辨率像元时也会出错。

为了避免这种不利因素的影响，因此，本发明进一步提出了采用第二增量△ndvis来对第一增量进行修正的构思。

其中，所述第二增量△ndvis指的是，采用现有envi/idl软件中集成的薄板样条(thinplatespline，tps)方法工具，直接将t1时间与tp时间的低空间分辨率modis像元分别插值到高空间分辨率的landsat像元上去(该插值方法工具已经集成在envi/idl软件中，为公知的方法，调用函数即可完成)。然后，使用tp时间的插值结果与t1时间的插值结果相减，其结果即为每个像元的所述第二增量△ndvis。由于该第二增量△ndvis忽略了地物类型的差异，不但可以提供空间上的信息，而且不受地物类型的干扰，因此通过下面的步骤将第一和第二增量进行融合，即可消除第一增量△ndvid缺陷和不利影响。

即，图1所示的心形标记的1个modis像元，对应下的64个landsat像元，每个像元的第一增量△ndvid与地物类型相关，图中表示有三种地物类型，则第1类地物类型下的24个像元的第一增量△ndvid都是相同的，第2类地物类型下的20个像元的第一增量△ndvid也是相同的，第3类地物类型下的20个像元的第一增量△ndvid也都是相同的，当然，这三种第一增量△ndvid的数值是分别计算获得的，是相互独立的。上述第一增量△ndvid是在假定同一像元的地物类别没有发生改变的情况下获得的，这显然会存在一定的问题。

而第二增量△ndvis的获取是根据薄板样条的方法工具直接将t1时间与tp时间的低空间分辨率的modis像元插值得到tp时间的高空间分辨率的landsat像元(薄板样条的方法原理是，在两张低空间分辨率modis图像中找出多个匹配点，应用薄板样条函数可以将这多个点形变到对应位置，同时给出了整个图像的插值)，将插值获得的tp时间的高空间分辨率的landsat图像的每个像元(例如对应于图1中的64个landsat像元)的ndvi数值与已知的t1时间的landsat图像的每个像元的ndvi数值相减，即可获得每个像元从t1至tp时间的第二增量△ndvis，即第二增量△ndvis与地物类型无关，每个像元均具有一个独立的第二增量△ndvis。当然，由于可以从低空间分辨率的modis像元直接插值得到高空间分辨率的landsat像元，则单独计算获得第二增量△ndvis是没有价值的，因此，此处提出第二增量△ndvis的目的就是为了下一步的融合，以消除单独采用第一增量△ndvid获得高时空分辨率ndvi时所带来的问题。

步骤b，使用贝叶斯模型平均方法，将每个高空间分辨率的像元在ndvi上的第一增量与第二增量进行权重组合计算，获得每个像元的ndvi综合增量。

这里使用两个权重来对每个像元(xi，yj)所对应的第一增量△ndvid(xi，yj)和第二增量△ndvis(xi，yj)进行权重组合。得到每个高空间分辨率像元(xi，yj)的ndvi综合增量为，

δndvicom(xi,yj)＝ws*δndvis(xi,yj)+wd*δndvid(xi,yj)(2)

公式(2)中，wd和ws分别表示第一增量和第二增量的权重，这两个权重是不随像元变化的，整幅图像只有一个wd和一个ws。△ndvicom(xi，yj)表示第一增量与第二增量结合得到的ndvi综合增量。

下面将需要使用贝叶斯模型平均(bayesianmodelaveraging，bma)方法来求解每个像元的这两个权重。bma方法是一种比较经典的多模型综合方法，它通过对似然函数进行优化来计算最终各个模型的权重。详细内容可以参看参考文献raftery，a.e.，gneiting，t.，balabdaoui，f.，&polakowski，m.(2005)，“使用贝叶斯模型平均方法对天气的集合预报校正”，每月天气评论，133，1155-1174(raftery，a.e.，gneiting，t.，balabdaoui，f.，&polakowski，m.(2005).usingbayesianmodelaveragingtocalibrateforecastensembles.monthlyweatherreview，133，1155-1174)。对于本发明中第一增量与第二增量的极大似然函数表示如下，

公式(3)中，l表示极大似然值，ws和wd分别表示两种增量的权重，△ndvid,i表示第i个像元的第一增量，△ndvis,i表示第i个像元的第二增量，△ndvii表示第i个像元的真实增量。gs和gd分别表示两个增量集合的概率密度函数，认为其为正态分布，n表示训练样本的数目。本发明中，由于没有高空间分辨率的真实ndvi增量数据，因此使用t1时间与tp时间低分辨率像元的真实增量作为训练样本，将高空间分辨率的第一增量和第二增量都进行8*8的合成升尺度到低空间分辨率上。这样即可得到低空间分辨率尺度的像元真实ndvi增量、低空间分辨率的第一增量和低空间分辨率的第二增量。根据低空间分辨的第一增量计算出其与低空间分辨率真实ndvi增量之间的平均平方误差(meansquareerror，mse)，这个平均平方误差就是概率密度函数gd的方差，同理可以获得概率密度函数gs的方差。gs(△ndvii|△ndvis,i)表示概率密度函数gs以其平均平方误差为方差，以△ndvis,i为均值，能产生△ndvii的概率。同样，gd(△ndvii|△ndvid,i)表示概率密度函数gd以其平均平方误差为方差，以△ndvid,i为均值，能产生△ndvii的概率。通过对似然函数(3)进行极大化，可以得到权重ws和wd。

值的说明的是，虽然ndvi综合增量△ndvicom(xi，yj)已经能够比较好的接近高空间分辨率像元从t1时间到tp时间的真实ndvi增量，但是其仍然可能包含误差，因此为了得到更准确的结果，在一个优选实施例中，还可以对获得的ndvi综合增量进行进一步的误差修正，其修正步骤为：

首先计算出低空间分辨率像元(x，y)整体的误差为，

公式(4)中，m表示低空间分辨率像元(x，y)所覆盖的高空间分辨率像元的数目，在图1中显示的即为64；r(x，y)表示低空间分辨率像元(x，y)整体的误差；△ndvi^c(x，y)表示tp时间与t1时间低空间分辨率像元(x，y)的ndvi的差值。

此处，需要将低空间分辨率像元尺度所包含的误差r(x，y)分配到其所覆盖的高空间分辨率像元(xi，yj)上。一般来说，影像上均质性比较高的区域(比如一大片水体或者一大片草地)往往预测时误差较小，而均质性低的区域(比如树木和裸土交错分布的稀疏林地)往往预测误差较大。因此，本发明根据高空间分辨率像元的均质性(xi，yj)来对误差r(x，y)进行分配。

即，然后计算高空间分辨率像元的均质性如下，

hi(xi,yj)＝nsmae/m，(5)

公式(5)中，hi(xi，yj)是以高空间分辨率像元(xi，yj)为中心，构建的8*8的窗口的均质性参数。m表示这个窗口中高空间分辨像元的数目，这里即64个，nsame表示这个窗口中，与中心像元(xi，yj)具有相同类型的高空间分辨率像元的数目。

之后，根据像元的均质性来对误差进行分配如下，

公式(6)中，r(xi，yj)即为高空间分辨率像元(xi，yj)应该分配到的误差。

最后，将这个误差加到公式(2)中的综合增量上去，得到修正过的高空间分辨率像元的ndvi增量如下，

公式(7)中，△ndvicom^new(xi，yj)即为最后计算得到的高空间分辨率像元(xi，yj)从t1时间到tp时间的修正过的ndvi综合增量。

步骤c，计算tp时间的高空间分辨率ndvi；

前面的步骤中已经计算得到的t1时间到tp时间，每个高空间分辨率像元(xi，yj)的ndvi综合增量，下面将综合增量加到t1时间的高空间分辨率ndvi上，即可得到tp时间的高空间分辨率ndvi数据。计算如下，

公式(8)中，ndvip(xi，yj)即为tp时间的高空间分辨率ndvi，ndvi1(xi，yj)即为t1时间的高空间分辨率ndvi。

如上所述，经过步骤a、b和c的处理后，预测时间tp所有的高空间分辨率像元都能通过时空加权的预测过程得到比较精确的融合结果。相比于已有的数据融合方法而言，本发明通过bma(贝叶斯模型平均)方法更充分的使用了低空间分辨率数据中所包含的信息，能够更好地应对土地覆盖类型变化造成的问题。

为了更好地说明本发明的技术效果，将本发明中的方法与背景技术部分提及的三种类型的数据融合方法进行比较，即针对starfm(空间和时间自适应反射融合模型)、ndvi-lmgm(ndvi-线性混合生长模型)和fsdaf(灵活时空数据融合方法)，分别进行时序数据的测试和单景数据的测试。测试中以ndvi为融合对象，低空间分辨率数据为modis数据，高空间分辨率数据为landsat数据。由于这两种数据获取的过境时间相同，并且轨道参数接近，因此比较适合于用来进行融合处理。

图2a-2f分别显示的是土地覆盖变化区域以及对应的landsatndvi数据和modisndvi数据；其中表示的某集水流域区的情况，具体选择的测试数据位于澳大利亚新南威尔士北部的gwydir下游流域研究区(lowergwydircatchmentsite)，经纬度分别为29°07′s和149°04′e。该区域是一个典型的农作物种植区，并于2004年12月发生了一次比较大的洪水事件，有大量农田被淹没，发生了比较明显的土地覆盖变化。为了验证本发明对于土地覆盖变化的检测效果，从澳大利亚联邦科学与工业研究组织(commonwealthscientificandindustrialresearchorganization，csiro)的官方网站：https://data.csiro.au/dap/landingpage？execution＝e2s2&_eventid＝viewdescription，下载到洪水发生前2004年11月26日和洪水发生后2004年12月12日的landsattm影像，对其进行了裁剪，计算出对应的ndvi(图2a和2b)。根据2004年11月26日裁剪过的影像(图2c)进行基于isodata的非监督分类，得到分类结果(图2f)。为了避免不同传感器之间的系统误差，以及几何纠正和大气校正等过程中带来的误差，本发明使用landsat数据进行8*8的升尺度合成得到modis数据(图2d和2e)。

图3a-3e分别显示的是几种不同的数据融合方法在土地覆盖变化区域的融合结果，以2004年11月26日的landsatndvi数据为参考数据，对2004年12月12日的数据进行融合。将四种方法的融合结果(图3a、图3b、图3c和图3d)与2004年12月12日真实的ndvi(图3e)进行比较。从局部放大图的整体可以看到，ndvi-lmgm方法的融合结果有明显的斑块问题，而本发明的方法融合的结果的平滑性比较好。从图3a-3e中较小的方框可以看到；starfm融合对于土地覆盖变化的检测能力不够，较小的方框中的黑色区域没有被识别出来，而本发明的方法能够很好地对土地覆盖的变化进行识别；从较大的方框可以看到，fsdaf融合的结果受到分类的影响相对大一些，融合结果中对于洪水发生前的水体边界有明显的保留，使得被水覆盖区域的过渡不平滑，而本发明的方法受到分类的影响相对小，对土地覆盖的识别能力更强，并且融合结果与真实值更接近。

综上所述，本发明基于贝叶斯模型平均方法将混合像元分解的增量与空间插值的增量进行综合，对增量的计算过程进行了优化。对时空数据进行了综合，有效提高了数据融合的精度。本发明充分利用了地理数据之间的时空关联性，为生产高时空分辨率的植被指数提供了一种行之有效的方法。

本领域技术人员应当理解，虽然本发明是按照多个实施例的方式进行描述的，但是并非每个实施例仅包含一个独立的技术方案。说明书中如此叙述仅仅是为了清楚起见，本领域技术人员应当将说明书作为一个整体加以理解，并将各实施例中所涉及的技术方案看作是可以相互组合成不同实施例的方式来理解本发明的保护范围。

以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的等同变化、修改与结合，均应属于本发明保护的范围。