*基于時空融合的NDVI時序生成技術在冬小麥監(jiān)測中的應用

李勝林，李大成，，韓啟金，龍小祥

(1a.太原理工大學 礦業(yè)工程學院，太原 030024；1b.太原理工大學 高分辨率對地觀測系統(tǒng)山西數(shù)據(jù)與應用中心，山西 晉中 030600；2.中國資源衛(wèi)星應用中心，北京 100094)

植被指數(shù)是根據(jù)植被在各波段的光譜反射特性進行一系列線性與非線性組合后得到的、生物生長過程中具有指示意義的數(shù)值[1]。植被指數(shù)可以從多個角度反映作物的生長環(huán)境和生長狀況，是利用遙感手段識別作物、監(jiān)測長勢和估計產(chǎn)量的重要參數(shù)。歸一化植被指數(shù)(normalized difference vegetation index，NDVI)是最常用的一種植被指數(shù)[2-3]，獲得高時序高空間分辨率的NDVI數(shù)據(jù)對冬小麥的生長監(jiān)測具有重要意義。

通過遙感數(shù)據(jù)可獲得高時序的NDVI數(shù)據(jù)，遙感數(shù)據(jù)可以在多時空尺度與空間范圍進行應用。然而，由于技術的限制以及像元的尺寸與掃描帶寬問題，目前尚沒有單一的傳感器能夠獲得既具有高空間分辨率又具有高時間分辨率的遙感數(shù)據(jù)[4]。高空間分辨率數(shù)據(jù)往往覆蓋范圍較小且重訪周期較長[5]，不適用于對植被的長期與實時監(jiān)測；而高時間分辨率遙感數(shù)據(jù)適用于大范圍實時監(jiān)測，但空間分辨率往往較低[6]。為了同時提高遙感數(shù)據(jù)的空間分辨率與時間分辨率，學者們提出了時空融合模型。

時空融合模型中影響最大的是由GAO et al[4]提出來的STARFM(spatial and temporal adaptive reflectance fusion model)算法。該算法利用同一時期的MODIS與Landsat反射率影像，結(jié)合預測時期的MODIS影像對預測時期的Landsat反射率影像進行預測，進而估算高時空分辨率的地表反射率；其預測精度較高，得到了廣泛的研究與應用[7-12]。ZHU et al[13]在STARFM的基礎上提出了ESTARFM(enhanced STARFM)，有效提高了在破碎地區(qū)STARFM的準確性；但是ESTARFM需要覆蓋相同研究區(qū)的兩景Landsat影像，不利于數(shù)據(jù)的獲取與大區(qū)域的應用。HUANG et al[14]最早利用稀疏表達與字典學習的方法提出SPSTFM(sparse-representation-based spatiotemporal reflectance fusion model)時空融合模型，用非線性的方式得到了包含高-低分辨率空間細節(jié)的字典對，使用字典對重建預測時刻的高分辨率影像；由于該模型計算復雜、運行效率較低，目前較少應用于大區(qū)域的研究。

STARFM時空融合模型并沒有限制數(shù)據(jù)源必須是Landsat與MODIS數(shù)據(jù)，但目前的相關研究中較少采用其他數(shù)據(jù)源，且算法參數(shù)考慮了國外衛(wèi)星數(shù)據(jù)特點，很少應用于國產(chǎn)衛(wèi)星數(shù)據(jù)。孫銳等[15]在STARFM的基礎上考慮了光譜差異，結(jié)合環(huán)境一號衛(wèi)星(HJ-1 CCD)數(shù)據(jù)與MODIS數(shù)據(jù)生成了研究區(qū)內(nèi)4-6月份的NDVI數(shù)據(jù)。MENG et al[16]于2011年在STARFM的基礎上提出了不同時空分辨率NDVI數(shù)據(jù)時空融合算法(spatial and temporal adaptive vegetation index fusion model,STAVFM)；于2013年利用STAVFM算法生成高時空NDVI數(shù)據(jù)，應用于作物的估產(chǎn)[17]。徐磊等[18]利用STAVFM算法，結(jié)合GF-1/WFV與MODIS數(shù)據(jù)實現(xiàn)了森林覆蓋定量提取。隨著我國高分一號衛(wèi)星(Gaofen-1 satellite,GF-1)的成功發(fā)射，利用高分一號衛(wèi)星16 m分辨率的多光譜寬覆蓋(Gaofen-1 satellite/wide field of view，GF-1/WFV)數(shù)據(jù)實現(xiàn)冬小麥生長動態(tài)監(jiān)測成為可能。理論上GF-1/WFV數(shù)據(jù)覆蓋周期為4 d，但實際中要達到4d的覆蓋周期需對某一地區(qū)進行高頻動態(tài)監(jiān)測，這樣勢必會大大增加GF-1衛(wèi)星的運行成本；再加上陰雨、云污染等天氣因素，僅通過GF-1/WFV數(shù)據(jù)獲得高時序的農(nóng)作物NDVI數(shù)據(jù)將會非常困難。采用時空融合算法將GF-1/WFV與MODIS數(shù)據(jù)進行融合得到高時序的GF-1/WFV數(shù)據(jù)成為解決此問題的重要途徑。

本文選取國產(chǎn)衛(wèi)星GF-1/WFV數(shù)據(jù)與MODIS反射率產(chǎn)品MOD09Q1數(shù)據(jù)，運用最為廣泛成熟的STARFM算法生成研究區(qū)內(nèi)冬小麥不同物候時期的GF-1/WFV NDVI數(shù)據(jù)，并對融合結(jié)果進行分析。采用JAKUBAUSKAS et al[19]提出的時間序列諧波分析法(harmonic analysis of time series,HANTS)進行濾波，以便減小噪聲影響，實現(xiàn)對GF-1/WFV NDVI時間序列數(shù)據(jù)更好的重構(gòu)。

1 研究區(qū)概況與數(shù)據(jù)預處理

1.1 研究區(qū)概況

研究區(qū)位于河南省東北部，新鄉(xiāng)市與鶴壁市交界處，研究區(qū)中心坐標為114°12′11″E，35°29′32″N，面積為19 km×19 km.如圖1所示，實驗區(qū)地處黃河下游，氣候溫和，地勢平坦，年平均氣溫13～15 ℃，年日照2 300～2 600 h，無霜期210～220 d，年均降水量600～800 mm，60%的降水量集中在6-8月份，春季易干旱，夏季易雨澇。實驗區(qū)主要地物為農(nóng)田(主要作物為小麥，還有少量玉米、棉花等)、居民住宅地以及少量裸地與道路等。

圖1 研究區(qū)位置示意圖Fig.1 Location of the study area

1.2 數(shù)據(jù)預處理

本文所用實驗數(shù)據(jù)為GF-1/WFV數(shù)據(jù)與MOD09Q1數(shù)據(jù)。獲取2015年11月初至2016年5月底研究區(qū)內(nèi)無云或者較少云量的GF-1/WFV數(shù)據(jù)15景，MOD09Q1數(shù)據(jù)31景。在進行融合處理前需要對GF-1/WFV數(shù)據(jù)與MOD09Q1數(shù)據(jù)進行預處理。對GF-1/WFV數(shù)據(jù)依次進行輻射定標、幾何校正、大氣校正，得到反射率數(shù)據(jù)，用于NDVI計算。輻射定標系數(shù)采用中國資源應用衛(wèi)星中心提供的2015年與2016年的定標系數(shù)。幾何校正以同一地區(qū)的Landsat數(shù)據(jù)為基準，以同名地物點為基礎進行二次多項式擬合。大氣校正使用ENVI自帶的FLASH模塊進行處理。對MOD09Q1數(shù)據(jù)采用MRT(MODIS reprojection tool)進行重新投影，得到與GF-1/WFV數(shù)據(jù)相同的坐標系統(tǒng)，并使用雙線性內(nèi)插法重采樣到16 m分辨率，然后進行幾何校正，以保證MOD09Q1與GF-1/WFV數(shù)據(jù)空間上匹配。

2 研究方法

2.1 STARFM算法

(1)

G(xi,yi,tk)=M(xi,yi,tk)+εk.

(2)

式中：G(xi,yi,tk)、M(xi,yi,tk)分別是GF-1/WFV和MOD09Q1數(shù)據(jù)在tk時的反射率值；εk是兩者的反射率差異(由帶寬和幾何畸變引起)。所需預測的t0時的G(xi,yi,t0)可表達如下：

G(xi,yi,t0)=M(xi,yi,t0)+ε0.

(3)

此算法進一步假設像元(xi,yi)處的地表覆蓋類型和系統(tǒng)誤差ε在t0和tk時不變，即ε0=εk.由此得出：

G(xi,yi,t0)=M(xi,yi,t0)+
G(xi,yi,tk)-M(xi,yi,tk) .

(4)

但實際上，上述理想情況并不能得到滿足，原因如下：MOD09Q1數(shù)據(jù)的掃描視場很大，其觀測量很可能并不是同質(zhì)像元；同一區(qū)域的地表覆蓋類型在兩個不同時間可能會發(fā)生改變；地表覆蓋狀態(tài)和由光照幾何變化所引起的雙向反射分布函數(shù)(BRDF)均有可能發(fā)生改變。

因此，STARFM算法利用一個權重函數(shù)Wijk來解決由上述問題所引起的反射率的預測誤差，用權重函數(shù)Wijk來計算t0時中心像元的反射率，最終得到預測時GF-1/WFV數(shù)據(jù)的反射率，即

(5)

其中，w是搜索框的尺寸，(xw/2,yw/2)是移動搜索框的中心像元。

STARFM算法根據(jù)光譜、時相和空間三個要素來計算最終的權重函數(shù)Wijk，主要包括光譜差異Sijk、時相差異Tijk、中心像元(xw/2,yw/2)和候選像元(xi,yi)在tk時的空間距離dijk.三要素表達式如下：

Sijk=|G(xi,yi,tk)-M(xi,yi,tk)|；

(6)

Tijk=|M(xi,yi,tk)-M(xi,yi,t0)|；

(7)

(8)

式中，Sijk值越小，細分辨率像元的光譜特征與周圍像元的差異就越小，該候選像元應分配較大的權重。Tijk越小說明在tk與t0之間植被變化越小，該候選像元應分配較大的權重。A是特定常數(shù)，dijk指的是實際的歐式距離，dijk值小的對應于較大的權重。聯(lián)合光譜差異、時相差異和空間距離可得到Cijk：

Cijk=dijk×Sijk×Tijk.

(9)

最后的權重可表示為：

(10)

本文利用該模型，集合GF-1/WFV數(shù)據(jù)與MOD09Q1數(shù)據(jù)，探究了STARFM算法對冬小麥關鍵物候期以及整個生長周期的監(jiān)測能力。

2.2 融合策略

本文基于MOD09Q1數(shù)據(jù)得到2015年10月中旬至2016年6月初研究區(qū)內(nèi)步長為8 d的MOD09Q1 NDVI數(shù)據(jù)(共31景)，并使用時間序列諧波分析法(HANTS)對研究區(qū)內(nèi)NDVI變化曲線進行平滑重建，如圖2所示。圖2中平滑后的時序NDVI曲線保持了平滑前NDVI曲線的大致形狀，能夠有效改善異常值，反映了研究區(qū)內(nèi)冬小麥生長期內(nèi)的動態(tài)變化與物候特征。

圖2 平滑前、后的MOD09Q1 NDVI曲線Fig.2 Unsmoothed and smoothed MOD09Q1 NDVI time series curves

通過分析圖2，可得到小麥生長的物候期。2015年10月24日左右小麥開始出苗，NDVI逐漸增大，至2015年12月19日左右NDVI值達到此階段峰值；之后進入越冬期，NDVI值緩慢下降。2016年2月2日左右，天氣回暖，小麥開始進入返青-拔節(jié)期，至2016年4月23日左右NDVI達到最大值。然后是抽穗到成熟的生殖生長階段，營養(yǎng)生長基本停止，一直到6月初小麥成熟，并伴隨著部分小麥的收割，在此期間NDVI值迅速下降。

本文以冬小麥不同物候期為研究對象，融合2015年11月1日與12月4日影像對冬小麥出苗生長期進行預測分析，融合2016年1月15日與2月8日影像對小麥越冬期進行預測分析，融合2016年3月11日與4月13日影像對小麥返青-拔節(jié)期進行預測分析，融合2016年4月13日與5月16日影像對小麥抽穗期進行預測分析，融合2016年5月4日與5月29日影像對小麥成熟期進行預測分析。鑒于一個月內(nèi)至少可以找到1景有效GF-1/WFV數(shù)據(jù)的事實，在時相上選擇基準日期與預測日期相差一個月左右進行預測分析，以便體現(xiàn)該算法的適用性。為考察時間間隔對預測精度的影響，分別以2016年1月15日、2月8日、3月11日、4月13日為基準影像，對剩下三個時相進行預測分析。研究技術路線圖如圖3所示。

圖3 研究技術路線圖Fig.3 Route chart of research technology

3 實驗結(jié)果與分析

圖4給出了不同物候期MOD09Q1與GF-1/WFV數(shù)據(jù)融合所得到的預測GF-1/WFV NDVI與該時相上實際GF-1/WFV NDVI散點圖的對比。使用平均絕對差值AAD(average absolute difference)、均方根誤差RMSE(root mean square error)、結(jié)構(gòu)相似性SSIM(structural similarity)、差分影像均值Mean和差分影像標準差SD(standard deviation)等評價指標進行分析，預測精度見表1.

由圖4可知，出苗生長期與成熟期預測精度較差，相關系數(shù)R分別為0.695 9、0.632 9，RMSE分別為0.100 7、0.168 4.越冬期、返青-拔節(jié)期、抽穗期具有較好的擬合精度：絕大部分散點比較一致地分布在直線y=x兩側(cè)；相關系數(shù)R均在0.840 0以上，表明預測NDVI與實際NDVI具有較好的空間一致性；RMSE均在0.066 3以下，離散程度較小，融合效果較好。同時對照表1，預測日期2016年2月8日、4月13日、5月16日的平均絕對差值(AAD<0.049 2)與結(jié)構(gòu)相似性(SSIM>0.858 0)比預測日期2015年12月4日(AAD=0.074 2，SSIM=0.712 7)與2016年5月29日(AAD=0.134 9，SSIM=0.505 4)有更好的預測結(jié)果。分析Mean與SD的變化情況，可以得到與上述一致的結(jié)論。

產(chǎn)生該結(jié)果的可能原因分析如下。以2015年11月1日預測2015年12月4日為例，這個時間段為小麥出苗與生長的快速階段，地表覆蓋類項可能會發(fā)生較大變化，因此地表覆蓋類項的變化嚴重影響了預測精度。以2016年5月4日預測2016年5月29日為例，這個時間段部分冬小麥已經(jīng)被收割，因此預測精度較差。以2016年1月15日預測2016年2月8日為例，此階段處于越冬期，小麥NDVI變化較小，地表覆蓋類型比較穩(wěn)定，但預測精度與以2016年3月11日預測2016年4月13日(此時間段內(nèi)NDVI變化較大)的預測精度比較一致，這說明NDVI的變化對預測精度并未產(chǎn)生較大影響。

表2列出了不同基準影像的預測結(jié)果精度分析，所選時間段內(nèi)地表類型比較穩(wěn)定。從表2可以看出：

物候期基準日期預測日期RRMSEAADSSIMSDMean出苗生長期2015-11-012015-12-040.695 90.100 70.074 20.712 70.120 20.118 0越冬期2016-01-152016-02-080.840 40.063 80.049 10.858 00.058 00.026 5返青-拔節(jié)期2016-03-112016-04-130.892 10.066 10.048 20.907 60.071 80.013 6抽穗期2016-04-132016-05-160.897 00.066 30.049 20.897 60.070 6-0.027 1成熟期2016-05-042016-05-290.632 90.168 40.034 90.505 40.118 00.120 2

表2 不同基準影像的預測精度Table 2 Prediction accuracy of different base image

1) 對于同一基準影像，隨著時間間隔的增加，各項評價指標顯示出預測精度趨于降低。這可能是因為隨著時間的增加，預測的不確定性增加得更多。

2) 預測精度因不同的預測方向而有一定差異。例如，以2016年1月15日預測2016年4月13日的結(jié)果與以2016年4月13日預測2016年1月15日的結(jié)果具有較大差異。

3) 地表類型穩(wěn)定時，物候引起的NDVI變化對預測結(jié)果影響較小。例如，以 2016年1月15日預測2016年4月13日，即使基準日期與預測日期相差3個月，相關系數(shù)也高至0.784 8，融合效果較好。

圖5為實際NDVI影像與預測NDVI影像的對比，可從中得到與圖4一致的結(jié)果。因此，在對冬小麥監(jiān)測時應充分考慮地表類型變化所帶來的差異，而物候變化引起的反射率變化對預測結(jié)果的影響較小。

左列圖為實際NDVI影像，右列圖為對應的預測NDVI影像圖5 實際GF-1/WFV NDVI影像與預測GF-1/WFV NDVI影像Fig.5 Observed NDVI and predicted NDVI from GF-1/WFV images

以MOD09Q1 8 d合成數(shù)據(jù)與研究區(qū)內(nèi)15景GF-1/WFV數(shù)據(jù)作為數(shù)據(jù)源，選擇與GF-1/WFV數(shù)據(jù)時相上最為接近的MOD09Q1數(shù)據(jù)作為基準影像，預測時選擇距預測日期最近的影像對作為基準影像對，得到研究區(qū)內(nèi)2015年10月8日至2016年6月2日間步長為8 d的GF-1/WFV NDVI數(shù)據(jù)，實現(xiàn)了對小麥生長周期內(nèi)的動態(tài)監(jiān)測。圖6給出了研究區(qū)內(nèi)實際MOD09Q1 NDVI均值與預測GF-1/WFV NDVI均值(均經(jīng)過HANTS濾波處理)，可看出二者具有較高的一致性。

總體而言，當基準影像與預測影像時間間隔在1個月以內(nèi)時，STARFM算法通過對MOD09Q1與GF-1/WFV的融合可以較好地擬合小麥的變化特征；但當?shù)乇砀采w類型發(fā)生變化時，預測精度較差。因此在冬小麥出苗時期與成熟期，應該增加更多時序數(shù)據(jù)以減小基準日期與預測日期的時間間隔，提高預測精度。

圖6 實際MOD09Q1 NDVI與預測GF-1/WFV NDVIFig.6 Observed MOD09Q1 NDVI and predicted GF-1/WFV NDVI

4 結(jié)束語

本文采用STARFM算法將MOD09Q1數(shù)據(jù)的時間變化信息與GF-1/WFV的空間差異進行結(jié)合，得到時間序列GF-1/WFV NDVI數(shù)據(jù)，有效地補充了缺失的GF-1/WFV NDVI數(shù)據(jù)，也證明了SATRFM算法在冬小麥監(jiān)測中具有一定應用潛力。在使用該算法生成高時序GF-1/WFV NDVI數(shù)據(jù)時，還需充分考慮到以下幾個影響因素：

1) 地表類型變化對預測結(jié)果帶來較大影響。在小麥出苗期與成熟收割期應增加影像對，減小基準日期與預測日期時間間隔，以得到更好的預測效果。

2) 當?shù)乇砀采w類型比較穩(wěn)定時，物候變化以及小麥生長對預測結(jié)果產(chǎn)生較小的影響；但隨著基準時間與預測時間間隔的增加，會有更多的不確定性因素，預測精度隨之降低。因此，預測時應盡量選擇較小時間間隔的影像對。

3) GF-1/WFV數(shù)據(jù)雖經(jīng)過大氣校正，但大氣條件及二向反射仍對其有一定影響，最終這些噪聲被引入到GF-1/WFV NDVI數(shù)據(jù)中。在以后的工作中應使用一種有效的去除噪聲機制，提高數(shù)據(jù)源質(zhì)量，減小預測的不確定性。

4) 不同傳感器之間的差異會造成影像的光譜差異。因此，以后應考慮到光譜響應函數(shù)，提取有效參數(shù)，減小不同傳感器影像間的光譜差異。