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

背景技術(shù)：

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

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

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

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

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

(3)基于混合模型的方法以fsdaf(靈活時(shí)空數(shù)據(jù)融合方法，flexiblespatiotemporaldatafusionmethod)為代表，參見zhu，x.，helmer，e.h.，gao，f.，liu，d.，chen，j.，&lefsky，m.a.(2016)，“一種靈活地將不同分辨率的衛(wèi)星影像進(jìn)行時(shí)空融合的方法”，環(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é)合起來，計(jì)算得到高分辨率影像上每個(gè)像元的增量。這類方法往往比較復(fù)雜，但是精度相對(duì)較高，而且能夠?qū)ν恋馗采w變化進(jìn)行一定程度的捕捉。

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

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

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

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

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

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

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

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

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

優(yōu)選地，所述步驟b中，每個(gè)高空間分辨率像元(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)種的每個(gè)像元的所述權(quán)重wd和ws通過貝葉斯模型平均方法求解獲得。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

附圖說明

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

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

圖2a-2f分別顯示的是土地覆蓋變化區(qū)域以及對(duì)應(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表示對(duì)應(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表示的是真實(shí)值。

具體實(shí)施方式

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

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

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

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

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

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

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

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

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

對(duì)應(yīng)于時(shí)間軸上方的1個(gè)modis像元的面積，等于時(shí)間軸下方的64個(gè)landsat像元的面積，將這64個(gè)landsat像元的ndvi增量按照其面積平均起來，即可認(rèn)定為等于時(shí)間軸上方的同樣面積的1個(gè)modis像元的ndvi增量。其計(jì)算過程為，根據(jù)高空間分辨率的landsatndvi數(shù)據(jù)的地物分類計(jì)算出每個(gè)低空間分辨率的modisndvi數(shù)據(jù)的像元中每種地物類型的豐度。分類數(shù)據(jù)可以來源于已有的分類產(chǎn)品，比如globeland30，也可以通過對(duì)高空間分辨率的landsatndvi數(shù)據(jù)的衛(wèi)星影像數(shù)據(jù)進(jìn)行分類得到。根據(jù)高空間分辨率的landsatndvi數(shù)據(jù)的地物分類結(jié)果，計(jì)算出每個(gè)低空間分辨率的modisndvi數(shù)據(jù)的像元中包含的不同地物類型的面積比例，即豐度。因此，將每種地物類型的ndvi增量與其豐度的乘積相加，就可以得到該面積下的全部landsat像元(64個(gè)像元)的ndvi增量的平均值，這個(gè)平均值正好等于前述同樣面積的modis像元(1個(gè)像元)從t1至tp時(shí)間的ndvi增量。

例如，圖1時(shí)間軸下方的64個(gè)landsat像元包括三種地物類型，則第1類地物類型的24個(gè)像元具有一個(gè)ndvi增量，第2類地物類型的20個(gè)像元具有另一個(gè)ndvi增量，第3類地物類型的20個(gè)像元具有又一個(gè)ndvi增量，這三種ndvi增量，相對(duì)于每個(gè)像元來說都是其從t1至tp時(shí)間的第一增量△ndvid。為了計(jì)算出這三個(gè)增量，需要構(gòu)建至少三個(gè)方程，此時(shí)假定最鄰近時(shí)間軸上方的心形標(biāo)記像元的8個(gè)像元(亦即圍繞心形標(biāo)記像元的8個(gè)像元)的ndvi增量與該心形標(biāo)記像元的ndvi增量相同，則可以以每個(gè)低空間分辨率像元為中心，構(gòu)建3*3的局部窗口。假設(shè)這個(gè)局部窗口中每一種地物類型在t1時(shí)間至tp時(shí)間內(nèi)的ndvi增量是相同的。根據(jù)這個(gè)局部窗口中的9個(gè)低空間分辨率像元的ndvi增量，可構(gòu)建如下的線性方程組，公式表示如下，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟c，計(jì)算tp時(shí)間的高空間分辨率ndvi；

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

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

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

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

圖2a-2f分別顯示的是土地覆蓋變化區(qū)域以及對(duì)應(yīng)的landsatndvi數(shù)據(jù)和modisndvi數(shù)據(jù)；其中表示的某集水流域區(qū)的情況，具體選擇的測(cè)試數(shù)據(jù)位于澳大利亞新南威爾士北部的gwydir下游流域研究區(qū)(lowergwydircatchmentsite)，經(jīng)緯度分別為29°07′s和149°04′e。該區(qū)域是一個(gè)典型的農(nóng)作物種植區(qū)，并于2004年12月發(fā)生了一次比較大的洪水事件，有大量農(nóng)田被淹沒，發(fā)生了比較明顯的土地覆蓋變化。為了驗(yàn)證本發(fā)明對(duì)于土地覆蓋變化的檢測(cè)效果，從澳大利亞聯(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影像，對(duì)其進(jìn)行了裁剪，計(jì)算出對(duì)應(yīng)的ndvi(圖2a和2b)。根據(jù)2004年11月26日裁剪過的影像(圖2c)進(jìn)行基于isodata的非監(jiān)督分類，得到分類結(jié)果(圖2f)。為了避免不同傳感器之間的系統(tǒng)誤差，以及幾何糾正和大氣校正等過程中帶來的誤差，本發(fā)明使用landsat數(shù)據(jù)進(jìn)行8*8的升尺度合成得到modis數(shù)據(jù)(圖2d和2e)。

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

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

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

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