采用地表反射率模型的Landsat時序數(shù)據集重構

李凱， 陳蕓芝， 汪小欽， 陳雪嬌

(1. 福州大學 空間數(shù)據挖掘與信息共享教育部重點實驗室， 福建 福州 350116；2. 福州大學 衛(wèi)星空間信息技術綜合應用國家地方聯(lián)合工程研究中心， 福建 福州 350116)

高時空分辨率遙感數(shù)據在監(jiān)測土地快速變化、作物生長,以及物候參數(shù)反演等方面具有重要的作用[1].Landsat系列衛(wèi)星以較長的數(shù)據、較高的空間分辨率，以及較好的傳感器性能在植被覆蓋變化監(jiān)測、植被生長評估等方面得到廣泛應用[2].但在多云地區(qū)，Landsat系列衛(wèi)星影像云覆蓋率較高，且晴空條件下獲取的影像數(shù)據時間間隔較長、季相不一致，這些都限制了影像的有效利用，使得這類傳感器難以實現(xiàn)高頻次、周期性的監(jiān)測.通過重構Landsat時序數(shù)據集，插補缺失觀測值，可加強對地表動態(tài)監(jiān)測的能力.噪聲的影響普遍存在于遙感圖像中，通過一系列方法降低時序數(shù)據中的噪聲水平，可為各種研究提供可靠的時序數(shù)據集[3].朱慧等[4]利用不同方法對重慶2010-2014年中分辨率成像光譜儀(MODIS)時序數(shù)據進行重建，認為WS(whittaker smoother)方法能較好地對不同土地覆蓋類型進行重建.宋春橋等[5]對藏北地區(qū)各植被類型擬合效果較好的非對稱高斯函數(shù)進行擬合，重構2007-2009年MODIS NDVI(歸一化植被指數(shù))時序數(shù)據集.劉亞南等[6]在傳統(tǒng)單Logistic模型的基礎上，采用參數(shù)構建模型對秦嶺樣區(qū)2001-2013年MODIS NDVI時序數(shù)據進行重建.李明等[7]利用同一類地物高質量像元均值代替噪聲的S-G濾波算法，對江西省2001-2003年MODIS NDVI時間序列進行重建.殷悅等[8]通過6種方法對鄱陽湖平原地區(qū)2001-2013年的地球觀測系統(tǒng)(SPOT)和MODIS時序數(shù)據進行重構，認為對于高植被覆蓋區(qū)域S-G濾波算法的去噪效果最優(yōu).為滿足大范圍、高精度、快速變化的遙感監(jiān)測，鄔明權等[9]提出多源遙感數(shù)據時空融合技術，此類方法可以克服中等空間分辨率影像時間分辨率過低的局限[10].郭文靜等[11]通過光譜差異改進ESTARFM(enhanced spatial and temporal adaptive reflectance fusion model)模型，對TM(thematic mapper)和AVHRR的NDVI數(shù)據進行融合，有效地實現(xiàn)若爾蓋高原NDVI數(shù)據集的重構.孫銳等[12]采用分類擬合的方法，計算光譜距離權重改進STARFM(spatial and temporal adaptive reflectance fusion model)模型，融合中國環(huán)境一號A/B星CCD(HJ-1CCD)，重構山東省西北部NDVI時序數(shù)據集.張猛等[13]利用STARFM算法，對Landsat OLI(operational land imager)和MODIS數(shù)據融合構建空間分辨率為30 m、時間間隔為16 d的時序數(shù)據集，用于區(qū)域植被凈初級生產力(NPP)的遙感估算.Meng等[14]基于STARFM算法，提出一種新的植被指數(shù)融合模型用于融合HJ-1CCD和MODIS數(shù)據，生成NDVI時序數(shù)據集.

目前，一系列基于降噪的時序數(shù)據重構方法已經取得了很好的成果，重構的遙感數(shù)據能準確地反映地表特征的變化[15]，但更適用于時間分辨率較高的時序數(shù)據集.時空融合算法利用時間信息和空間信息進行時序數(shù)據集的重構，但存在一定的局限性.首先，時空融合技術非常依賴較低空間分辨率影像的觀測值；其次，該方法始終需要清晰無云的影像.因此，對于快速發(fā)生變化的區(qū)域，使用時空融合的方法進行變化監(jiān)測效率較低.針對目前時間序列數(shù)據重構存在的局限性，Zhu等[16]提出地表反射率模型，該模型不依賴高時間分辨率影像，生成任意時間的清晰無云的Landsat影像.本文以福州市區(qū)某東部面積近300 km2的區(qū)域作為研究區(qū)，采用地表反射率模型進行Landsat時序數(shù)據集重構.

1 研究區(qū)概況與數(shù)據源

1.1 研究區(qū)概況

研究區(qū)位于福建省福州市區(qū)東部，北緯為26°0′55.81″～26°8′47.41″，東經為119°21′45.16″～119°33′39.22″，面積約為300 km2.研究區(qū)位置，如圖1所示.

圖1 研究區(qū)位置Fig.1 Position of study area

由圖1可知：地形以山地丘陵為主，植被類型多樣，森林覆蓋率達到57.8%.研究區(qū)屬典型的亞熱帶季風氣候，全年云雨天氣較多，質量好的遙感影像相對較少.

1.2 數(shù)據源及其預處理

基于研究目標，選取2013-2016年云覆蓋90%以下的Landsat OLI影像，收集到的影像共44景.通過對原始影像進行大氣校正，消除大氣吸收和散射造成的輻射量誤差.大部分影像受云、云陰影等影響嚴重，對于少量云覆蓋的影像，采用Landsat OLI自帶的第9波段云檢測，設定閾值對云進行掩膜；對于含云量較高的影像，采用基于對象的Fmask算法進行云和云陰影的檢測.由于去云、去云陰影后影像依然存在異常值，擬合去云、去云陰影后的時序影像.將擬合影像與Fmask算法檢測結果進行對比，結果表明：對于沒有檢測出的異常值，其在藍光波段的值明顯高于模型的擬合值[17].

2 研究方法

時序地表反射率預測模型由諧波模型和一個趨勢分量組成，用于預測每個波段的地表反射率，模型較為穩(wěn)定、強健，不容易受噪聲(云，云陰影等)影響.Zhu等[16,18]利用地表反射率模型,重構美國新英格蘭地區(qū)1982-2011年Landsat影像，并利用擬合影像進行土地覆蓋的連續(xù)變化監(jiān)測；利用1982-2013年Landsat TM和ETM+影像，重構美國6個地區(qū)的Landsat時序數(shù)據集，并對重構影像進行評價.

地表反射率模型由4個系數(shù)組成:a0,i用于估算i波段整體的像元值或平均值；a1,i,b1,i用于模擬i波段由物候和太陽角度差異引起的年內變化；c1,i表示i波段的年際變化，用于估計i波段的長期趨勢.影像去云、去云陰影和去異常值等預處理后，剩余每個像元的有效觀測值個數(shù)達到12個，保證模型的精度和穩(wěn)定.對于有效觀測值個數(shù)小于12大于6的像元，采用地表反射率模型；對于有效觀測值小于5個的像元，使用所有有效觀測值的中位數(shù)來表示地表反射率.

采用最小二乘法(OLS)擬合，時序模型為

上式中:x為日期；T為一年的天數(shù)(T=365.25).

3 結果與分析

3.1 擬合結果分析

為研究區(qū)內每個Landsat像元生成時間序列模型，重構2013-2016年時序數(shù)據集，對2013-2016年所有可用的Landsat OLI觀測值擬合NDVI時間序列.選取荷葉、闊葉樹、柳樹、馬尾松4種不同類型的植被像元，擬合結果與原始NDVI對比，如圖2所示.圖2中：時間序列編號1～16代表2013年，17～39代表2014年，40～61代表2015年，62～84代表2016年.

由圖2可知：原始序列中大部分點的NDVI值都與擬合序列較為接近；對于季節(jié)特征變化大的荷葉，NDVI值變化范圍明顯大于其他3種植被，對于四季常綠的闊葉樹和馬尾松NDVI值全年保持在較高的水平；隨著有效像元數(shù)的增加，模型擬合效果越好，擬合結果與原始NDVI值越為接近.因此，擬合結果與原始有效像元保持較好的一致性，模型能夠擬合不同植被穩(wěn)定的生長變化過程.

(a) 荷葉 (b) 闊葉樹

(c) 柳樹 (d) 馬尾松圖2 擬合結果與原始NDVI對比Fig.2 Comparison of fitting results with original NDVI

3.2 實測光譜驗證

2016年7月22日、2016年11月16日和2016年12月23日分別對研究區(qū)內荷葉、闊葉樹、柳樹、馬尾松4種類型植被的光譜(圖1)進行測量，利用實測數(shù)據，驗證模型擬合結果.

擬合結果與實測NDVI對比，如圖3所示.由圖3可知：闊葉樹、柳樹的擬合結果變化趨勢與實測NDVI基本保持一致，7-11月份NDVI值逐漸增大，11-12月份緩慢降低；馬尾松、荷葉的NDVI值變化幅度小，NDVI保持在較高的水平，7-11月份的擬合結果變化幅度小于實測NDVI變化幅度，11-12月份兩者變化趨勢基本保持一致；樣點的實測值高于擬合結果的NDVI值，主要是因為在像元尺度下，存在混合信息現(xiàn)象，導致擬合結果的NDVI值偏低，而實測值為一個純凈點且植被本身存在覆蓋不均勻.因此，不同類型植被的擬合結果與實測值保持了較好的一致性.

(a) 荷葉 (b) 闊葉樹

(c) 柳樹 (d) 馬尾松 圖3 擬合結果與實測NDVI對比Fig.3 Comparison of fitting results with measured NDVI

3.3 原始影像驗證

選擇不同季節(jié)且影像質量較好的原始影像與擬合影像進行均方根誤差分析，擬合結果與原始影像均方根誤差，如表1所示.表1中：RMSE為均方根誤差；Blue，Green，Red，NIR為影像波段的名稱.

由表1可知：可見光波段的RMSE較小(以地表反射率為單位，約為0.01);近紅外波段的RMSE值較大，這是因為研究區(qū)的大部分被植被覆蓋，植被在近紅外波段能產生較高的反射率;擬合結果與原始影像非常相近，差異很小，擬合精度較高.

分析2015年第270天的擬合結果與原始影像各波段相關性，擬合結果與原始影像散點圖，如圖4所示.圖4中：x為擬合結果；y為原始影像;R為相關系數(shù).由圖4可知：擬合結果與原始影像的散點分布集中，4個波段都集中在直線y=x附近；各波段的相關系數(shù)均高于0.9，表明兩者具有較強的相關性,整體擬合精度較高.

(a) Blue (b) Green

(c) Red (d) NIR圖4 擬合結果與原始影像散點圖Fig.4 Fitting results and original image scatter image

為了進一步驗證模型的能力，擬合4個不同季節(jié)(2015年第270天；2016年第49天；2016年第97天；2016年第209天)影像，擬合結果與原始影像對比圖，如圖5所示.由圖5可知：Landsat OLI波段5，4，3合成標準假彩色圖像，植被在影像中大致呈紅色；2015年第270天與2016年第209天的擬合影像紅色區(qū)域面積明顯高于2016年第49天與2016年第97天，夏秋兩季植被覆蓋明顯高于春冬兩季;擬合影像空間信息與地表覆蓋的變化過程均與同一時間的原始影像保持較高的一致性，并且能夠清晰地表達植被的物候差異；對于真實Landsat影像被云和云陰影覆蓋的地方，擬合影像仍然可以提供清晰的觀測.

(a) 2015年第270天(原始)(b) 2016年第49天(原始)(c) 2016年第97天(原始)(d) 2016年第209天(原始)

(e) 2015年第270天(擬合)(f) 2016年第49天(擬合)(g) 2016年第97天(擬合)(h) 2016年第209天(擬合)圖5 擬合結果與原始影像對比圖Fig.5 Comparison of fitting results with original image

擬合結果與原始影像NDVI差值圖，如圖6所示.由圖6可知：2015年第270天、2016年第49天、2016年第209天這3個時間點的大部分區(qū)域NDVI差異較小，極小部分差異較大的區(qū)域主要集中在地表覆蓋變化幅度大的區(qū)域；圖5(c)影像存在云的影響，原始影像紅光波段的反射率整體偏高，導致擬合結果與原始影像NDVI的差值偏大；擬合結果與原始影像的差異較小，說明模型的擬合效果較好.

(a) 2015年第270天 (b) 2016年第49天

(c) 2016年第97天 (d) 2016年第209天圖6 擬合結果與原始影像差值圖Fig.6 Difference of fitting results with original image

4 結論

采用地表反射率模型，基于2013-2016年云覆蓋90%以下的所有Landsat OLI影像，構建Landsat遙感時序數(shù)據集.通過地表反射率模型獲得的擬合結果與真實影像之間整體具有較高的相關性，各波段的相關系數(shù)均高于0.9，均方根誤差較低，可見光波段在0.01左右，近紅外波段在0.02左右.對不同類型植被均能取得較好的擬合效果，不同類型植被擬合NDVI與實測NDVI變化趨勢基本保持一致.擬合結果能清晰表現(xiàn)出植被的物候差異，與原始影像之間的差異較小，其空間信息及地表覆蓋的變化過程也與原始影像保持較高的一致性.通過地表反射率模型提供的密集清晰無云影像可用于植被動態(tài)變化監(jiān)測研究，進一步提高植被變化監(jiān)測的準確性和時效性.

該模型需要大量清晰的Landsat OLI影像進行精確的時間序列模型預測，對于長期多云和雪覆蓋的區(qū)域，可能沒有足夠的清晰影像滿足模型的要求，無法進行精確的時間序列模型預測，下一步將針對這方面進行改進，提出更加全面的時序重建方法.