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

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

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

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

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

1.1 研究區(qū)概況

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

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

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

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

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

2 研究方法

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

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

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

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

3 結(jié)果與分析

3.1 擬合結(jié)果分析

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

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

(a) 荷葉 (b) 闊葉樹

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

3.2 實(shí)測(cè)光譜驗(yàn)證

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

擬合結(jié)果與實(shí)測(cè)NDVI對(duì)比，如圖3所示.由圖3可知：闊葉樹、柳樹的擬合結(jié)果變化趨勢(shì)與實(shí)測(cè)NDVI基本保持一致，7-11月份NDVI值逐漸增大，11-12月份緩慢降低；馬尾松、荷葉的NDVI值變化幅度小，NDVI保持在較高的水平，7-11月份的擬合結(jié)果變化幅度小于實(shí)測(cè)NDVI變化幅度，11-12月份兩者變化趨勢(shì)基本保持一致；樣點(diǎn)的實(shí)測(cè)值高于擬合結(jié)果的NDVI值，主要是因?yàn)樵谙裨叨认拢嬖诨旌闲畔F(xiàn)象，導(dǎo)致擬合結(jié)果的NDVI值偏低，而實(shí)測(cè)值為一個(gè)純凈點(diǎn)且植被本身存在覆蓋不均勻.因此，不同類型植被的擬合結(jié)果與實(shí)測(cè)值保持了較好的一致性.

(a) 荷葉 (b) 闊葉樹

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

3.3 原始影像驗(yàn)證

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

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

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

(a) Blue (b) Green

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

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

(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 擬合結(jié)果與原始影像對(duì)比圖Fig.5 Comparison of fitting results with original image

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

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

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

4 結(jié)論

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

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