本發(fā)明涉及一種多源遙感數(shù)據(jù)融合的方法，主要用來生產(chǎn)高時空分辨率的歸一化植被指數(shù)(ndvi)，屬于時空數(shù)據(jù)融合的范疇。

背景技術(shù)：

歸一化植被指數(shù)(normalizeddifferencevegetationindex，縮寫為ndvi)是一種廣泛使用的植被指數(shù)。ndvi數(shù)值能夠反映植被的生長狀況并與植被的生理參數(shù)比如葉綠素含量，含水量以及葉面積指數(shù)等有直接的聯(lián)系。

cn102831310b公開了一種構(gòu)建高空間分辨率的時間序列數(shù)據(jù)的方法，其背景技術(shù)部分對于ndvi數(shù)據(jù)的相關(guān)知識進行了詳細說明，其中提及兩種現(xiàn)有的ndvi數(shù)據(jù)：landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)。其中，landsatndvi數(shù)據(jù)中的ndvi數(shù)值的空間分辨率為30m*30m，時間分辨率為16天，而modisndvi數(shù)據(jù)中的ndvi數(shù)值的最大空間分辨率為250m*250m，時間分辨率為1天。也就是說，landsatndvi數(shù)據(jù)具有較高的空間分辨率和較低的時間分辨率，而modisndvi數(shù)據(jù)具有較低的空間分辨率和較高的時間分辨率。

為了能夠同時獲得高空間分辨率和高時間分辨率的ndvi數(shù)據(jù)，現(xiàn)有技術(shù)發(fā)展出了多種對landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)進行融合，以獲取高空間分辨率的ndvi數(shù)據(jù)的方法，主要包括三種類型：(1)基于權(quán)函數(shù)的方法；(2)基于混合像元分解的方法；(3)基于混合模型的方法。

其中，(1)基于權(quán)函數(shù)的方法以starfm(空間和時間自適應(yīng)反射融合模型，spatialandtemporaladaptivereflectancefusionmodel)為典型代表，參見gao，f.，masek，j.，schwaller，m.，&hall，f.(2006)，“對landsat與modis的地表反射率進行融合：預(yù)測每日的landsat地表反射率”，電氣與電子工程師學(xué)會之地球科學(xué)與遙感學(xué)報，44，2207-2218。(gao，f.，masek，j.，schwaller，m.，&hall，f.(2006).ontheblendingofthelandsatandmodissurfacereflectance:predictingdailylandsatsurfacereflectance.ieeetransactionsongeoscienceandremotesensing，44，2207-2218)。該類方法著眼于相似像元的搜索及其權(quán)重的確定，通過當(dāng)前像元與周圍像元在時間、空間以及光譜等特征上的相似性確定相似像元并計算對應(yīng)的權(quán)重，以相似像元的平均增量作為當(dāng)前像元的增量。這類方法往往在純凈度比較高的區(qū)域才有比較好的效果。

(2)基于混合像元分解的方法以ndvi-lmgm(ndvi-線性混合生長模型，ndvilinearmixinggrowthmodel)為代表，參見rao，y.，zhu，x.，chen，j.，&wang，j.(2015)，“一種升級過地使用多時相modisndvi數(shù)據(jù)與landsattm/etm+影像生產(chǎn)高時空分辨率的ndvi時序數(shù)據(jù)的方法”，遙感，7，7865-7891(rao，y.，zhu，x.，chen，j.，&wang，j.(2015).animprovedmethodforproducinghighspatial-resolutionndvitimeseriesdatasetswithmulti-temporalmodisndvidataandlandsattm/etm+images.remotesensing，7，7865-7891)。該類方法通過混合像元分解，計算出每個像元在時間上的增量，從而得到預(yù)測時間的ndvi數(shù)據(jù)，有關(guān)該方法的具體內(nèi)容也可參見前述現(xiàn)有技術(shù)cn102831310b，二者實質(zhì)上是相同的。這類方法對于異質(zhì)性程度比較高的區(qū)域也能夠有比較好的效果，但其假設(shè)土地覆蓋狀況不發(fā)生變化，這一點不一定能夠滿足。

(3)基于混合模型的方法以fsdaf(靈活時空數(shù)據(jù)融合方法，flexiblespatiotemporaldatafusionmethod)為代表，參見zhu，x.，helmer，e.h.，gao，f.，liu，d.，chen，j.，&lefsky，m.a.(2016)，“一種靈活地將不同分辨率的衛(wèi)星影像進行時空融合的方法”，環(huán)境遙感，172，165-177(zhu，x.，helmer，e.h.，gao，f.，liu，d.，chen，j.，&lefsky，m.a.(2016).aflexiblespatiotemporalmethodforfusingsatelliteimageswithdifferentresolutions.remotesensingofenvironment，172，165-177)。該類方法將多種模型，如權(quán)函數(shù)模型、混合像元分解模型以及空間插值模型等結(jié)合起來，計算得到高分辨率影像上每個像元的增量。這類方法往往比較復(fù)雜，但是精度相對較高，而且能夠?qū)ν恋馗采w變化進行一定程度的捕捉。

以上三類方法在各自適用的領(lǐng)域里都取得了比較好的效果。但是它們都沒有考慮到使用兩期低空間分辨率數(shù)據(jù)插值到高空間分辨率的方式，兩期插值結(jié)果之間的增量能夠從一定程度上反映像元尺度的增量，同時也包含了土地覆蓋變化的信息。本發(fā)明基于貝葉斯模型平均方法將混合像元分解的增量與空間插值的增量進行綜合，能夠有效提高數(shù)據(jù)融合的精度。為生產(chǎn)高時空分辨率的ndvi提供了新的途徑。

技術(shù)實現(xiàn)要素：

本發(fā)明要解決的技術(shù)問題是提供一種基于時空加權(quán)的生產(chǎn)高時空分辨率ndvi的方法，以減少或避免前面所提到的問題。

為解決上述技術(shù)問題，本發(fā)明提出了一種基于時空加權(quán)的生產(chǎn)高時空分辨率ndvi的方法，用于將高空間分辨率和低時間分辨率的landsatndvi數(shù)據(jù)與低空間分辨率和高時間分辨率的modisndvi數(shù)據(jù)融合，從基準時間t1獲得預(yù)測時間tp的高時空分辨率的ndvi數(shù)據(jù)，所述方法包括如下步驟：

步驟a：以同一基準時間t1為起點，利用modisndvi數(shù)據(jù)和landsatndvi數(shù)據(jù)，分別獲得預(yù)測時間tp的每個像元在ndvi上的第一增量△ndvid和第二增量△ndvis；其中，

所述第一增量△ndvid指的是，從基準時間t1到預(yù)測時間tp，假定modisndvi數(shù)據(jù)的每個像元的ndvi增量，等于同樣面積下的landsatndvi數(shù)據(jù)的不同地物類別的像元的ndvi增量的平均值，則利用t1時間的modisndvi數(shù)據(jù)和landsatndvi數(shù)據(jù)，以及tp時間的modisndvi數(shù)據(jù)，計算獲得tp時間到t1時間的landsatndvi數(shù)據(jù)的每個像元的ndvi增量，該增量即為所述第一增量△ndvid；

所述第二增量△ndvis指的是，采用現(xiàn)有envi/idl軟件中集成的薄板樣條方法工具，直接將t1時間與tp時間的低空間分辨率modis像元分別插值到高空間分辨率的landsat像元上去，然后使用tp時間的插值結(jié)果與t1時間的插值結(jié)果相減，其結(jié)果即為每個像元的所述第二增量△ndvis；

步驟b，使用貝葉斯模型平均方法，將步驟a獲得的每個高空間分辨率的像元在ndvi上的所述第一增量與所述第二增量進行權(quán)重組合計算，獲得每個像元的ndvi綜合增量；

步驟c，通過步驟b計算得到的t1時間到tp時間，每個像元的所述ndvi綜合增量，加到t1時間對應(yīng)的高空間分辨率ndvi數(shù)值上，即可得到tp時間的高空間分辨率ndvi數(shù)值。

優(yōu)選地，所述步驟b中，每個高空間分辨率像元(xi，yj)的所述ndvi綜合增量用公式表示為，

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

其中，公式(2)中，wd和ws分別表示所述第一增量和所述第二增量的權(quán)重，△ndvicom(xi，yj)表示所述第一增量與所述第二增量結(jié)合得到的所述ndvi綜合增量。

優(yōu)選地，所述公式(2)種的每個像元的所述權(quán)重wd和ws通過貝葉斯模型平均方法求解獲得。

優(yōu)選地，所述步驟b進一步包括對所述ndvi綜合增量進行誤差修正的步驟。

優(yōu)選地，所述誤差修正的步驟包括：

首先計算出低空間分辨率像元(x，y)整體的誤差為，

公式(4)中，m表示低空間分辨率像元(x，y)所覆蓋的高空間分辨率像元的數(shù)目；r(x，y)表示低空間分辨率像元(x，y)整體的誤差；△ndvi^c(x，y)表示tp時間與t1時間低空間分辨率像元(x，y)的ndvi的差值；

然后計算高空間分辨率像元的均質(zhì)性如下，

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

公式(5)中，hi(xi，yj)是以高空間分辨率像元(xi，yj)為中心，構(gòu)建的8*8的窗口的均質(zhì)性參數(shù)；m表示這個窗口中高空間分辨像元的數(shù)目，nsame表示這個窗口中，與中心像元(xi，yj)具有相同類型的高空間分辨率像元的數(shù)目；

之后，根據(jù)像元的均質(zhì)性來對誤差進行分配如下，

公式(6)中，r(xi，yj)即為高空間分辨率像元(xi，yj)應(yīng)該分配到的誤差；

最后，將這個誤差加到公式(2)中的綜合增量上去，得到修正過的高空間分辨率像元的ndvi增量如下，

公式(7)中，△ndvicom^new(xi，yj)即為最后計算得到的高空間分辨率像元(xi，yj)從t1時間到tp時間的修正過的所述ndvi綜合增量。

優(yōu)選地，所述步驟c中，利用所述修正過的ndvi綜合增量計算tp時間的高空間分辨率ndvi的公式為：

公式(8)中，ndvip(xi，yj)即為tp時間的每個像元(xi，yj)的高空間分辨率ndvi，ndvi1(xi，yj)即為t1時間的每個像元(xi，yj)的高空間分辨率ndvi。

本發(fā)明通過對低空間分辨率的數(shù)據(jù)和高空間分辨率數(shù)據(jù)進行混合像元分解以及空間插值得到高空間分辨率數(shù)據(jù)在時間上的兩種增量，然后采用貝葉斯模型平均方法將所述兩種增量進行權(quán)重組合計算，得到優(yōu)化過的綜合增量，提高了增量計算的準確性，充分利用了地理數(shù)據(jù)間的時空聯(lián)系，提高了融合精度。

附圖說明

以下附圖僅旨在于對本發(fā)明做示意性說明和解釋，并不限定本發(fā)明的范圍。其中，

圖1顯示的是根據(jù)本發(fā)明的一個具體實施例的第一增量的原理示意圖；

圖2a-2f分別顯示的是土地覆蓋變化區(qū)域以及對應(yīng)的landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)；其中表示的某集水流域區(qū)的情況，分別為：圖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表示對應(yīng)的非監(jiān)督分類結(jié)果。

圖3a-3e分別顯示的是幾種不同的數(shù)據(jù)融合方法在土地覆蓋變化區(qū)域的融合結(jié)果；其中圖3a-3d分別表示的是2004年12月12日landsatndvi經(jīng)過ndvi-lmgm融合、starfm融合、fsdaf融合、本發(fā)明的方法融合的結(jié)果，圖3e表示的是真實值。

具體實施方式

為了對本發(fā)明的技術(shù)特征、目的和效果有更加清楚的理解，現(xiàn)對照附圖說明本發(fā)明的具體實施方式。但本領(lǐng)域的技術(shù)人員應(yīng)該知道，以下實施例并不是對本發(fā)明技術(shù)方案作的唯一限定，凡是在本發(fā)明技術(shù)方案精神實質(zhì)下所做的任何等同變換或改動，均應(yīng)視為屬于本發(fā)明的保護范圍。

下面詳細說明根據(jù)本發(fā)明提供的基于時空加權(quán)的生產(chǎn)高時空分辨率ndvi的方法的原理，當(dāng)然，與背景部分所述的三種類型的對landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)進行融合的方法的一樣，本發(fā)明也是對上述不同數(shù)據(jù)源的遙感數(shù)據(jù)進行融合的方法，亦即，本發(fā)明的基于時空加權(quán)的生產(chǎn)高時空分辨率ndvi的方法用于將高空間分辨率和低時間分辨率的landsatndvi數(shù)據(jù)與低空間分辨率和高時間分辨率的modisndvi數(shù)據(jù)融合，從基準時間t1獲得預(yù)測時間tp的高時空分辨率的ndvi數(shù)據(jù)。

由于獲取landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)的衛(wèi)星有相似的軌道參數(shù)以及不足30分鐘的衛(wèi)星過境時間間隔(該間隔對于每天24小時來說可以忽略不計)，因此，在理論上可以假設(shè)這兩種數(shù)據(jù)同一地點和同一基準時間t1的ndvi數(shù)據(jù)僅存在空間分辨率上的差別，其余均可忽略不計。其中，landsatndvi數(shù)據(jù)的空間分辨率為30m*30m，時間分辨率為16天；modisndvi數(shù)據(jù)的空間分辨率為250m*250m，時間分辨率為1天；則數(shù)據(jù)融合之后生產(chǎn)獲得的高時空分辨率的ndvi數(shù)據(jù)，其最大空間分辨率可以等于landsatndvi數(shù)據(jù)的空間分辨率(30m*30m)，其最大時間分辨率可以等于modisndvi數(shù)據(jù)的時間分辨率(1天)。也就是說，通過已知空間分辨率為250m*250m、時間分辨率為1天的modisndvi數(shù)據(jù)，以及空間分辨率為30m*30m、時間分辨率為16天的landsatndvi數(shù)據(jù)，以同一基準時間t1為起點，本發(fā)明的方法可以預(yù)測獲得其后15天的任意1天的預(yù)測時間tp的高時空分辨率的ndvi數(shù)據(jù)。

下面詳細說明本發(fā)明的基于時空加權(quán)的生產(chǎn)高時空分辨率ndvi的方法的具體步驟。

步驟a：以同一基準時間t1為起點，利用modisndvi數(shù)據(jù)和landsatndvi數(shù)據(jù)，分別獲得預(yù)測時間tp的每個像元在ndvi上的第一增量△ndvid和第二增量△ndvis。

其中，所述第一增量△ndvid指的是，從基準時間t1到預(yù)測時間tp，假定modisndvi數(shù)據(jù)的每個像元的ndvi增量，等于同樣面積下的landsatndvi數(shù)據(jù)的不同地物類別的像元的ndvi增量的平均值，則利用t1時間的modisndvi數(shù)據(jù)和landsatndvi數(shù)據(jù)，以及tp時間的modisndvi數(shù)據(jù)，計算獲得tp時間到t1時間的landsatndvi數(shù)據(jù)的每個像元的ndvi增量，該增量即為所述第一增量△ndvid。

具體來說，基準時間t1(實際上存在30分鐘的時間間隔)的低空間分辨率modisndvi數(shù)據(jù)和高空間分辨率的landsatndvi數(shù)據(jù)對同一地點來說僅存在空間分辨率上的差別，為了融合時將兩種數(shù)據(jù)的位置對準方便計算，可以將低空間分辨率的modisndvi數(shù)據(jù)從250m*250m重采樣到240*240m，從而一個modis像元可以完全覆蓋8*8＝64個landsat像元。

如圖1所示，其顯示的是根據(jù)本發(fā)明的一個具體實施例的第一增量的原理示意圖，其中時間軸的上方表示的是低空間分辨率的modisndvi數(shù)據(jù)，時間軸下方表示的是高空間分辨率的landsatndvi數(shù)據(jù)。以同一基準時間t1為起點，利用t1時間的modisndvi數(shù)據(jù)和landsatndvi數(shù)據(jù)，以及tp時間的modisndvi數(shù)據(jù)，可以獲得時間軸右下方虛線表示的原本不存在的tp時間的高空間分辨率的landsatndvi數(shù)據(jù)。

此時，假定從基準時間t1到預(yù)測時間tp，modisndvi數(shù)據(jù)的每個像元的ndvi增量，等于同樣面積下的landsatndvi數(shù)據(jù)的不同地物類別的像元的ndvi增量的平均值。即，圖1中時間軸上方心形表示的1個modis像元從t1至tp時間的ndvi增量，可以通過將tp時間的低空間分辨率的modisndvi數(shù)據(jù)與t1時間的低空間分辨率的modisndvi數(shù)據(jù)直接相減得到。

對應(yīng)于時間軸上方的1個modis像元的面積，等于時間軸下方的64個landsat像元的面積，將這64個landsat像元的ndvi增量按照其面積平均起來，即可認定為等于時間軸上方的同樣面積的1個modis像元的ndvi增量。其計算過程為，根據(jù)高空間分辨率的landsatndvi數(shù)據(jù)的地物分類計算出每個低空間分辨率的modisndvi數(shù)據(jù)的像元中每種地物類型的豐度。分類數(shù)據(jù)可以來源于已有的分類產(chǎn)品，比如globeland30，也可以通過對高空間分辨率的landsatndvi數(shù)據(jù)的衛(wèi)星影像數(shù)據(jù)進行分類得到。根據(jù)高空間分辨率的landsatndvi數(shù)據(jù)的地物分類結(jié)果，計算出每個低空間分辨率的modisndvi數(shù)據(jù)的像元中包含的不同地物類型的面積比例，即豐度。因此，將每種地物類型的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。為了計算出這三個增量，需要構(gòu)建至少三個方程，此時假定最鄰近時間軸上方的心形標記像元的8個像元(亦即圍繞心形標記像元的8個像元)的ndvi增量與該心形標記像元的ndvi增量相同，則可以以每個低空間分辨率像元為中心，構(gòu)建3*3的局部窗口。假設(shè)這個局部窗口中每一種地物類型在t1時間至tp時間內(nèi)的ndvi增量是相同的。根據(jù)這個局部窗口中的9個低空間分辨率像元的ndvi增量，可構(gòu)建如下的線性方程組，公式表示如下，

公式(1)中，△ndvi^c(x，y)表示從t1到tp時間低空間分辨率像元(x，y)的ndvi增量，fl(x，y)表示第l種類型在低空間分辨率像元(x，y)中的豐度，△ndvil(x，y)表示3*3的局部窗口中第l種類型的ndvi增量。

根據(jù)公式(1)通過最小二乘進行求解可以得到△ndvil(x，y)。對于圖示局部窗口中的低空間分辨率像元(x，y)而言，其覆蓋范圍下的高分辨率像元(xi，yj)在時間上的增量為△ndvid(xi，yj)，這個增量即稱為第一增量，下標(i，j)表示低空間分辨率像元覆蓋下的第(i，j)個高空間分辨率像元。已經(jīng)假設(shè)這個局部窗口中，每種地物類型的ndvi增量是固定的。因此，對于高空間分辨率像元(xi，yj)而言，其在時間上的增量△ndvid(xj，yj)就等于高空間分辨率像元(xj，yj)所屬的那個地物類型的ndvi增量。如果高空間分辨率像元(xi，yj)屬于第l類，那么△ndvid(xi，yj)就等于△ndvil(x，y)。

對于第一增量△ndvid而言，由于tp時間的地物類型與t1時間的地物類型可能有較大不同，例如水土流失、森林砍伐、火災(zāi)、洪水等均可造成不同時間的地物類型發(fā)生變化。因此，公式(1)中的每種地物類型的豐度fl(x，y)在兩個時間下有可能是不同的，會導(dǎo)致地物類型的ndvi增量△ndvil(x，y)計算出錯。而且在將每種地物類型的ndvi增量△ndvil(x，y)按類型分配給tp時間的高空間分辨率像元時也會出錯。

為了避免這種不利因素的影響，因此，本發(fā)明進一步提出了采用第二增量△ndvis來對第一增量進行修正的構(gòu)思。

其中，所述第二增量△ndvis指的是，采用現(xiàn)有envi/idl軟件中集成的薄板樣條(thinplatespline，tps)方法工具，直接將t1時間與tp時間的低空間分辨率modis像元分別插值到高空間分辨率的landsat像元上去(該插值方法工具已經(jīng)集成在envi/idl軟件中，為公知的方法，調(diào)用函數(shù)即可完成)。然后，使用tp時間的插值結(jié)果與t1時間的插值結(jié)果相減，其結(jié)果即為每個像元的所述第二增量△ndvis。由于該第二增量△ndvis忽略了地物類型的差異，不但可以提供空間上的信息，而且不受地物類型的干擾，因此通過下面的步驟將第一和第二增量進行融合，即可消除第一增量△ndvid缺陷和不利影響。

即，圖1所示的心形標記的1個modis像元，對應(yīng)下的64個landsat像元，每個像元的第一增量△ndvid與地物類型相關(guān)，圖中表示有三種地物類型，則第1類地物類型下的24個像元的第一增量△ndvid都是相同的，第2類地物類型下的20個像元的第一增量△ndvid也是相同的，第3類地物類型下的20個像元的第一增量△ndvid也都是相同的，當(dāng)然，這三種第一增量△ndvid的數(shù)值是分別計算獲得的，是相互獨立的。上述第一增量△ndvid是在假定同一像元的地物類別沒有發(fā)生改變的情況下獲得的，這顯然會存在一定的問題。

而第二增量△ndvis的獲取是根據(jù)薄板樣條的方法工具直接將t1時間與tp時間的低空間分辨率的modis像元插值得到tp時間的高空間分辨率的landsat像元(薄板樣條的方法原理是，在兩張低空間分辨率modis圖像中找出多個匹配點，應(yīng)用薄板樣條函數(shù)可以將這多個點形變到對應(yīng)位置，同時給出了整個圖像的插值)，將插值獲得的tp時間的高空間分辨率的landsat圖像的每個像元(例如對應(yīng)于圖1中的64個landsat像元)的ndvi數(shù)值與已知的t1時間的landsat圖像的每個像元的ndvi數(shù)值相減，即可獲得每個像元從t1至tp時間的第二增量△ndvis，即第二增量△ndvis與地物類型無關(guān)，每個像元均具有一個獨立的第二增量△ndvis。當(dāng)然，由于可以從低空間分辨率的modis像元直接插值得到高空間分辨率的landsat像元，則單獨計算獲得第二增量△ndvis是沒有價值的，因此，此處提出第二增量△ndvis的目的就是為了下一步的融合，以消除單獨采用第一增量△ndvid獲得高時空分辨率ndvi時所帶來的問題。

步驟b，使用貝葉斯模型平均方法，將每個高空間分辨率的像元在ndvi上的第一增量與第二增量進行權(quán)重組合計算，獲得每個像元的ndvi綜合增量。

這里使用兩個權(quán)重來對每個像元(xi，yj)所對應(yīng)的第一增量△ndvid(xi，yj)和第二增量△ndvis(xi，yj)進行權(quán)重組合。得到每個高空間分辨率像元(xi，yj)的ndvi綜合增量為，

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

公式(2)中，wd和ws分別表示第一增量和第二增量的權(quán)重，這兩個權(quán)重是不隨像元變化的，整幅圖像只有一個wd和一個ws?！鱪dvicom(xi，yj)表示第一增量與第二增量結(jié)合得到的ndvi綜合增量。

下面將需要使用貝葉斯模型平均(bayesianmodelaveraging，bma)方法來求解每個像元的這兩個權(quán)重。bma方法是一種比較經(jīng)典的多模型綜合方法，它通過對似然函數(shù)進行優(yōu)化來計算最終各個模型的權(quán)重。詳細內(nèi)容可以參看參考文獻raftery，a.e.，gneiting，t.，balabdaoui，f.，&polakowski，m.(2005)，“使用貝葉斯模型平均方法對天氣的集合預(yù)報校正”，每月天氣評論，133，1155-1174(raftery，a.e.，gneiting，t.，balabdaoui，f.，&polakowski，m.(2005).usingbayesianmodelaveragingtocalibrateforecastensembles.monthlyweatherreview，133，1155-1174)。對于本發(fā)明中第一增量與第二增量的極大似然函數(shù)表示如下，

公式(3)中，l表示極大似然值，ws和wd分別表示兩種增量的權(quán)重，△ndvid,i表示第i個像元的第一增量，△ndvis,i表示第i個像元的第二增量，△ndvii表示第i個像元的真實增量。gs和gd分別表示兩個增量集合的概率密度函數(shù)，認為其為正態(tài)分布，n表示訓(xùn)練樣本的數(shù)目。本發(fā)明中，由于沒有高空間分辨率的真實ndvi增量數(shù)據(jù)，因此使用t1時間與tp時間低分辨率像元的真實增量作為訓(xùn)練樣本，將高空間分辨率的第一增量和第二增量都進行8*8的合成升尺度到低空間分辨率上。這樣即可得到低空間分辨率尺度的像元真實ndvi增量、低空間分辨率的第一增量和低空間分辨率的第二增量。根據(jù)低空間分辨的第一增量計算出其與低空間分辨率真實ndvi增量之間的平均平方誤差(meansquareerror，mse)，這個平均平方誤差就是概率密度函數(shù)gd的方差，同理可以獲得概率密度函數(shù)gs的方差。gs(△ndvii|△ndvis,i)表示概率密度函數(shù)gs以其平均平方誤差為方差，以△ndvis,i為均值，能產(chǎn)生△ndvii的概率。同樣，gd(△ndvii|△ndvid,i)表示概率密度函數(shù)gd以其平均平方誤差為方差，以△ndvid,i為均值，能產(chǎn)生△ndvii的概率。通過對似然函數(shù)(3)進行極大化，可以得到權(quán)重ws和wd。

值的說明的是，雖然ndvi綜合增量△ndvicom(xi，yj)已經(jīng)能夠比較好的接近高空間分辨率像元從t1時間到tp時間的真實ndvi增量，但是其仍然可能包含誤差，因此為了得到更準確的結(jié)果，在一個優(yōu)選實施例中，還可以對獲得的ndvi綜合增量進行進一步的誤差修正，其修正步驟為：

首先計算出低空間分辨率像元(x，y)整體的誤差為，

公式(4)中，m表示低空間分辨率像元(x，y)所覆蓋的高空間分辨率像元的數(shù)目，在圖1中顯示的即為64；r(x，y)表示低空間分辨率像元(x，y)整體的誤差；△ndvi^c(x，y)表示tp時間與t1時間低空間分辨率像元(x，y)的ndvi的差值。

此處，需要將低空間分辨率像元尺度所包含的誤差r(x，y)分配到其所覆蓋的高空間分辨率像元(xi，yj)上。一般來說，影像上均質(zhì)性比較高的區(qū)域(比如一大片水體或者一大片草地)往往預(yù)測時誤差較小，而均質(zhì)性低的區(qū)域(比如樹木和裸土交錯分布的稀疏林地)往往預(yù)測誤差較大。因此，本發(fā)明根據(jù)高空間分辨率像元的均質(zhì)性(xi，yj)來對誤差r(x，y)進行分配。

即，然后計算高空間分辨率像元的均質(zhì)性如下，

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

公式(5)中，hi(xi，yj)是以高空間分辨率像元(xi，yj)為中心，構(gòu)建的8*8的窗口的均質(zhì)性參數(shù)。m表示這個窗口中高空間分辨像元的數(shù)目，這里即64個，nsame表示這個窗口中，與中心像元(xi，yj)具有相同類型的高空間分辨率像元的數(shù)目。

之后，根據(jù)像元的均質(zhì)性來對誤差進行分配如下，

公式(6)中，r(xi，yj)即為高空間分辨率像元(xi，yj)應(yīng)該分配到的誤差。

最后，將這個誤差加到公式(2)中的綜合增量上去，得到修正過的高空間分辨率像元的ndvi增量如下，

公式(7)中，△ndvicom^new(xi，yj)即為最后計算得到的高空間分辨率像元(xi，yj)從t1時間到tp時間的修正過的ndvi綜合增量。

步驟c，計算tp時間的高空間分辨率ndvi；

前面的步驟中已經(jīng)計算得到的t1時間到tp時間，每個高空間分辨率像元(xi，yj)的ndvi綜合增量，下面將綜合增量加到t1時間的高空間分辨率ndvi上，即可得到tp時間的高空間分辨率ndvi數(shù)據(jù)。計算如下，

公式(8)中，ndvip(xi，yj)即為tp時間的高空間分辨率ndvi，ndvi1(xi，yj)即為t1時間的高空間分辨率ndvi。

如上所述，經(jīng)過步驟a、b和c的處理后，預(yù)測時間tp所有的高空間分辨率像元都能通過時空加權(quán)的預(yù)測過程得到比較精確的融合結(jié)果。相比于已有的數(shù)據(jù)融合方法而言，本發(fā)明通過bma(貝葉斯模型平均)方法更充分的使用了低空間分辨率數(shù)據(jù)中所包含的信息，能夠更好地應(yīng)對土地覆蓋類型變化造成的問題。

為了更好地說明本發(fā)明的技術(shù)效果，將本發(fā)明中的方法與背景技術(shù)部分提及的三種類型的數(shù)據(jù)融合方法進行比較，即針對starfm(空間和時間自適應(yīng)反射融合模型)、ndvi-lmgm(ndvi-線性混合生長模型)和fsdaf(靈活時空數(shù)據(jù)融合方法)，分別進行時序數(shù)據(jù)的測試和單景數(shù)據(jù)的測試。測試中以ndvi為融合對象，低空間分辨率數(shù)據(jù)為modis數(shù)據(jù)，高空間分辨率數(shù)據(jù)為landsat數(shù)據(jù)。由于這兩種數(shù)據(jù)獲取的過境時間相同，并且軌道參數(shù)接近，因此比較適合于用來進行融合處理。

圖2a-2f分別顯示的是土地覆蓋變化區(qū)域以及對應(yīng)的landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)；其中表示的某集水流域區(qū)的情況，具體選擇的測試數(shù)據(jù)位于澳大利亞新南威爾士北部的gwydir下游流域研究區(qū)(lowergwydircatchmentsite)，經(jīng)緯度分別為29°07′s和149°04′e。該區(qū)域是一個典型的農(nóng)作物種植區(qū)，并于2004年12月發(fā)生了一次比較大的洪水事件，有大量農(nóng)田被淹沒，發(fā)生了比較明顯的土地覆蓋變化。為了驗證本發(fā)明對于土地覆蓋變化的檢測效果，從澳大利亞聯(lián)邦科學(xué)與工業(yè)研究組織(commonwealthscientificandindustrialresearchorganization，csiro)的官方網(wǎng)站：https://data.csiro.au/dap/landingpage？execution＝e2s2&_eventid＝viewdescription，下載到洪水發(fā)生前2004年11月26日和洪水發(fā)生后2004年12月12日的landsattm影像，對其進行了裁剪，計算出對應(yīng)的ndvi(圖2a和2b)。根據(jù)2004年11月26日裁剪過的影像(圖2c)進行基于isodata的非監(jiān)督分類，得到分類結(jié)果(圖2f)。為了避免不同傳感器之間的系統(tǒng)誤差，以及幾何糾正和大氣校正等過程中帶來的誤差，本發(fā)明使用landsat數(shù)據(jù)進行8*8的升尺度合成得到modis數(shù)據(jù)(圖2d和2e)。

圖3a-3e分別顯示的是幾種不同的數(shù)據(jù)融合方法在土地覆蓋變化區(qū)域的融合結(jié)果，以2004年11月26日的landsatndvi數(shù)據(jù)為參考數(shù)據(jù)，對2004年12月12日的數(shù)據(jù)進行融合。將四種方法的融合結(jié)果(圖3a、圖3b、圖3c和圖3d)與2004年12月12日真實的ndvi(圖3e)進行比較。從局部放大圖的整體可以看到，ndvi-lmgm方法的融合結(jié)果有明顯的斑塊問題，而本發(fā)明的方法融合的結(jié)果的平滑性比較好。從圖3a-3e中較小的方框可以看到；starfm融合對于土地覆蓋變化的檢測能力不夠，較小的方框中的黑色區(qū)域沒有被識別出來，而本發(fā)明的方法能夠很好地對土地覆蓋的變化進行識別；從較大的方框可以看到，fsdaf融合的結(jié)果受到分類的影響相對大一些，融合結(jié)果中對于洪水發(fā)生前的水體邊界有明顯的保留，使得被水覆蓋區(qū)域的過渡不平滑，而本發(fā)明的方法受到分類的影響相對小，對土地覆蓋的識別能力更強，并且融合結(jié)果與真實值更接近。

綜上所述，本發(fā)明基于貝葉斯模型平均方法將混合像元分解的增量與空間插值的增量進行綜合，對增量的計算過程進行了優(yōu)化。對時空數(shù)據(jù)進行了綜合，有效提高了數(shù)據(jù)融合的精度。本發(fā)明充分利用了地理數(shù)據(jù)之間的時空關(guān)聯(lián)性，為生產(chǎn)高時空分辨率的植被指數(shù)提供了一種行之有效的方法。

本領(lǐng)域技術(shù)人員應(yīng)當(dāng)理解，雖然本發(fā)明是按照多個實施例的方式進行描述的，但是并非每個實施例僅包含一個獨立的技術(shù)方案。說明書中如此敘述僅僅是為了清楚起見，本領(lǐng)域技術(shù)人員應(yīng)當(dāng)將說明書作為一個整體加以理解，并將各實施例中所涉及的技術(shù)方案看作是可以相互組合成不同實施例的方式來理解本發(fā)明的保護范圍。

以上所述僅為本發(fā)明示意性的具體實施方式，并非用以限定本發(fā)明的范圍。任何本領(lǐng)域的技術(shù)人員，在不脫離本發(fā)明的構(gòu)思和原則的前提下所作的等同變化、修改與結(jié)合，均應(yīng)屬于本發(fā)明保護的范圍。